Hi - i'm trying to model an alignment of four genes, but assuming a tree and branch lengths, with a separate GTR model for each partition. To do this, I want to constrain the overall mean rate to be equal to the branch lengths of the trees I give it.. This all seems to work fine, but I can't work out how to constrain the mean rate of my four partitions.
The rates of the four partitions are given by gene1rate, gene2rate etc., and I want to constrain the sum of rate*(numsite/totalsites) to be one.
I initially tried to do the following, but I think I realise what i'm doing, rather than constraining the sum of meanrates to be 1.0, the line overallrate := 1.0 clears the dependency of overallrate := meanrate1 + meanrate2 etc. and instead enforces that overallrate is 1! This seems obvious in retrospect.
How should I construct this slightly complicated constraint, given that the := operator only allows a single variable on the LHS? I'm sure there's an obvious way to do it, but i can't work it out.. I'm guessing you've come across similar problems and will see it straight away!
Thanks,
James
Code:DataSetFilter filteredData12S = CreateFilter (ds,1,DATA_FILE_PARTITION_MATRIX[1][0]);
DataSetFilter filteredDataPOS1 = CreateFilter (ds,1,DATA_FILE_PARTITION_MATRIX[1][1]);
DataSetFilter filteredDataPOS2 = CreateFilter (ds,1,DATA_FILE_PARTITION_MATRIX[1][2]);
DataSetFilter filteredDataPOS3 = CreateFilter (ds,1,DATA_FILE_PARTITION_MATRIX[1][3]);
ACCEPT_BRANCH_LENGTHS = 1;
ACCEPT_ROOTED_TREES = 1;
global gene1rate = 1;
global gene2rate = 1;
global gene3rate = 1;
global gene4rate = 1;
/*constrain all these to 1.*/
global totalsize := filteredData12S.sites + filteredDataPOS1.sites + filteredDataPOS2.sites + filteredDataPOS3.sites;
fprintf (stdout,"DATASET IS ",totalsize," CHARS\n");
fprintf (stdout,"PARTITIONS ARE ",filteredData12S.sites," ",filteredDataPOS1.sites," ",filteredDataPOS2.sites," ",filteredDataPOS3.sites,"\n");
global meanrate1 := gene1rate * (filteredData12S.sites / totalsize);
global meanrate2 := gene2rate * (filteredDataPOS1.sites / totalsize);
global meanrate3 := gene3rate * (filteredDataPOS2.sites / totalsize);
global meanrate4 := gene4rate * (filteredDataPOS3.sites / totalsize);
global overallrate := meanrate1 + meanrate2 + meanrate3 + meanrate4;
global overallrate := 1.0;