Hi,

Does HyPhy handle unconstrained, non-time reversible models? I'm trying to implement an 11 rate-parameter model (similar to PAML nucleotide model 8), but so far I've had not luck - the results seem very suspicious. Most likely I'm doing something drastically wrong.

I want the rate parameters to be constant across the phylogeny, so I define my parameters and matrix as below:

ACCEPT_ROOTED_TREES=1;

global a;
global b;
global c;
global d;
global e;
global f;
global g;
global h;
global i;
global j;
global k;
global l;

UNRESTRateMatrix = {{*,m*a,m*b,m*c}{m*d,*,m*e,m*f}{m*g,m*h,*,m*i}{m*j,m*l,m*k,*}};
Model UNREST = (UNRESTRateMatrix, obsFreqs, 0);

Here are some problems:

1) I don't want to multiply the rate matrix by nucleotide frequencies, since all the parameters can vary independently and I don't need time-reversibility - > I use the 0 option as the third parameter of  the "Model". Why do I still need to supply the frequency vector as the second argument? How are those frequencies used? The likelihood, tree, and parameter values change if I use equalFreqs vs. obsFreqs, even though I don't multiply the matrix by the frequency vector.

2. The model, as it is above, seems over-parametriized. I think I should be able to get rid of one of the global rate parameters (a - k). I tried this by setting one of the parameters to 1 (i.e. replacing, say, a by 1) in the matrix definition, but this alters the likelihood and is dependent on which parameter is held constant.

What am I doing wrong and how should I parametrize the rate matrix? Thanks in advance for any help,

Jacek

Hi Sergei.

Thanks for the reply. I like the second option - solving the equations analytically for the equilibrium frequencies, since that would allow me to use the unconstrained model and interpret the parameters directly as relative rates of substitutions. BUT, I remember at some point trying to solve the equilibrium distribution for the 12 parameter model, and the solution is painfully long...

So I chose the first suggestion: I used the NRM.mdl definition. However, I seem to run into the same problem that I had before. When I use your original code, which sets AG as 1, I get an answer that for my dataset may differ by as much as 20 LogLikelihood units from setting another parameter, GT = 1 (I made sure that I replaced GT = 1 in both the matrix definition and the global parameter definition for TG).

     global TA:=AT+(AC-CA)*r1__+(AG-GA)*r2__;
     global TC:=CT+(CA-AC)*r0__+(CG-GC)*r2__;
     global TG:=1+(GA-AG)*r0__+(GC-CG)*r1__;

     
     ModelMatrixName =
                 {{*,t*AC,t*AG,t*AT}
                  {t*CA,*,t*CG,t*CT}
                  {t*GA,t*GC,*,t*1}
                  {t*TA,t*TC,t*TG,*}};
                 
Fixing most of the other parameters alters the likelihood to some extent.

I also tried this with the GTR model, and I still get differences when using different parameters as a reference (although not as big, about 2 units).

I've redone this many times now, and can't find a problem. Do you have any ideas?

Jacek

Thanks!!!

I'm glad I asked instead of banging my head against the wall. The MAXIMUM_ITERATIONS_PER_VARIABLE was indeed the problem. Now everything converges fine.  I would very much appreciate  and example of how to:

"deduce correct stationary distribution for a model run-time (so you don't need to solve analytically)".

HyPhy seems like a great tool, but expectedly with all its flexibility it takes a while to figure out. I may be bothering you again...

Jacek

Hello,

In the code that is provided here, is the output the substitution rates? If not, how can I calculate them?

Thanks,

Marissa

Hello Sergei,

Thank you for the help and the prompt response!

After thinking about this more, I think I was not asking about quite the right thing. To give you an idea about what it is that I am looking for, I will explain what I am doing and ask if you have any insight about what I should be calculating.

I'm trying to compare the nucleotide substitution process in two different data sets. Both of these consist of ancestor descendent alignments where we have an ancestral transposon aligned to the human genome. I would like to be able to compare the rate of one particular substitution (say C->T) in the two data sets.

The processes I'm interested in are not time-reversible, so the model you created in your 2006 code earlier in this thread seems to be the right one to use.

I believe that what I want to compare is the instantaneous rate matrix for the two data sets. However, it seems to me that with the model in your 2006 code, I cannot compare the values I get for C to T substitutions (in the rate matrix) between two different data sets, because everything is relative to A to G substitutions. Is that correct?

Is the instantaneous rate matrix what I want? And if so, is there a way to compare it between the two datasets?

Thanks,

Marissa

Dear Sergei,

For our two data sets, we know that the evolutionary time is the same for the two data sets. Given that information, is there some way to compare the rates?

Also, here is an attempt I made at doing what you suggest (which to me sounds like calculating a joint likelihood and performing a log ratio test), modified from your 2006 code. However, I think this code isn't doing quite what I want it to. I would certainly appreciate any pointers!

[code]
/* Load data 1 */
DataSet seq1 = ReadDataFile ("Data/human/NNNA2.fa");
DataSetFilter filter1 = CreateFilter (seq1, 1);

/* Load data 2 */
DataSet seq2 = ReadDataFile ("Data/human/NNNG2.fa");
DataSetFilter filter2 = CreateFilter (seq2, 1);

/* Define general non-reversible model */
/* This file defines the transition matrix for the General Non-Reversible model
  without explicit frequency parameterization.
  01/04/2006  by Sergei L. Kosakovsky Pond
*/
global AC = 1;
global AG = 1;
global AT = 1;
global CA = 1;
global CG = 1;
global CT = 1;
global GA = 1;
global GC = 1;
global GT = 1;
global TA = 1;
global TC = 1;
global TG = 1;

function PopulateModelMatrix (ModelMatrixName&)
{
     ModelMatrixName =
                 {{*,t*AC,t,t*AT}
                  {t*CA,*,t*CG,t*CT}
                  {t*GA,t*GC,*,t*GT}
                  {t*TA,t*TC,t*TG,*}};
     return 0;
}

/* NRM = Non-Reversible Model */
NRM = 0;
MULTIPLY_BY_FREQS = PopulateModelMatrix ("NRM");

/* Calculate equilibrium frequencies */
dynamicFreqs = {4,1};
vectorOfFrequencies = {4,1};
vectorOfFrequencies [0] := computeEQF (AC,AT,CG,CT,GT,CA,GA,TA,GC,TC,TG);
vectorOfFrequencies [1] := dynamicFreqs[1];
vectorOfFrequencies [2] := dynamicFreqs[2];
vectorOfFrequencies [3] := dynamicFreqs[3];

Model GRMModel = (NRM, vectorOfFrequencies, MULTIPLY_BY_FREQS);
FREQUENCY_SENSITIVE = 0;
ACCEPT_ROOTED_TREES = 1;

/* Create two trees, one for seq1 and one for seq2.
  We set the anc branch to 0 by setting t to 0. */
Tree tree1 = (anc,der);
Tree tree2 = (anc,der);
tree1.anc.t := 0;
tree2.anc.t := 0;

/*---------- Unconstrained Model ----------*/

LikelihoodFunction theLnLik = (filter1, tree1, filter2, tree2);

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

/*----------- Constrained Model -----------*/

LikelihoodFunction theLnLik = (filter1, tree1, filter2, tree2);

/* Constraint */
tree1.der.CT := tree2.der.CT;

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


/* Compute and print the log ratio */
LRT = 2 (unconstrLik - constrainedLnLik);
fprintf  (stdout, "LRT: ", LRT, "\n");


/*--------------- function to compute equilibrium frequencies ---------------*/

ffunction computeEQF (dummy1,dummy2,dummy3,dummy4,dummy5,dummy6,dummy7,dummy8,dummy9,dummy10,dummy11)
{
     t = 1;
     tempMatrix = NRM;
     for (k=0; k<4; k=k+1)
     {
           tempMatrix[k][3] = 1;
     }
     tempMatrix = Inverse(tempMatrix);
     for (k=0; k<4; k=k+1)
     {
           dynamicFreqs[k] = tempMatrix[3][k];
     }
     return dynamicFreqs[0];
}
[/code]

Thank you for all the help!

Marissa

Dear Sergei,

I fixed that up and now it's working fine. Thanks!

Now that I know that the instantaneous substitution rate (in my case CT) is larger in one data set than the other, what I would also like to know is the magnitude of the difference between the C->T rates between the data sets, given we know that the evolutionary time is the same for both data sets.

In your 2006 code, it looks like I should be comparing the values in NRM (and in my fixed code, NRM_2 for the second data set). However, I did some testing with modeled data, and the values in NRM seem to fluctuate related to the ACGT content of the sequence even when the underlying substitution rates are the same, but not in any obvious way that I could simply factor out.

My question in, what do the values in NRM represent and what should I be comparing to find out the relative substitution rates between my two data sets?

Thanks in advance,

Marissa

Here is the code I've been using:

[code]
/* Load data 1 */
DataSet seq1 = ReadDataFile ("Data/modeled/AT1_s_1.fa");
DataSetFilter filter1 = CreateFilter (seq1, 1);

/* Load data 2 */
DataSet seq2 = ReadDataFile ("Data/human/AT2_s_1.fa");
DataSetFilter filter2 = CreateFilter (seq2, 1);

/* Define general non-reversible model */
/* This file defines the transition matrix for the General Non-Reversible model
  without explicit frequency parameterization.
  01/04/2006  by Sergei L. Kosakovsky Pond
*/

global AC = 1;
global AG = 1;
global AT = 1;
global CA = 1;
global CG = 1;
global CT = 1;
global GA = 1;
global GC = 1;
global GT = 1;
global TA = 1;
global TC = 1;
global TG = 1;

global AC_2 = 1;
global AG_2 = 1;
global AT_2 = 1;
global CA_2 = 1;
global CG_2 = 1;
global CT_2 = 1;
global GA_2 = 1;
global GC_2 = 1;
global GT_2 = 1;
global TA_2 = 1;
global TC_2 = 1;
global TG_2 = 1;


function PopulateModelMatrix (ModelMatrixName&)
{
     ModelMatrixName =
                 {{*,t*AC,t,t*AT}
                  {t*CA,*,t*CG,t*CT}
                  {t*GA,t*GC,*,t*GT}
                  {t*TA,t*TC,t*TG,*}};
     return 0;
}

function PopulateModelMatrix_2 (ModelMatrixName&)
{
     ModelMatrixName =
                 {{*, t_2*AC_2, t_2, t_2*AT_2}
                  {t_2*CA_2, *, t_2*CG_2, t_2*CT_2}
                  {t_2*GA_2, t_2*GC_2, *, t_2*GT_2}
                  {t_2*TA_2, t_2*TC_2, t_2*TG_2, *}};
     return 0;
}


/* Create two trees, one for seq1 and one for seq2.
  We set the anc branch to 0 by setting t to 0. */

/* NRM = Non-Reversible Model */
NRM = 0;
MULTIPLY_BY_FREQS = PopulateModelMatrix ("NRM");

/* Calculate equilibrium frequencies */
dynamicFreqs = {4,1};
vectorOfFrequencies = {4,1};
vectorOfFrequencies [0] := computeEQF (AC,AT,CG,CT,GT,CA,GA,TA,GC,TC,TG);
vectorOfFrequencies [1] := dynamicFreqs[1];
vectorOfFrequencies [2] := dynamicFreqs[2];
vectorOfFrequencies [3] := dynamicFreqs[3];

/* Create the model for the */
Model GRMModel = (NRM, vectorOfFrequencies, MULTIPLY_BY_FREQS);
FREQUENCY_SENSITIVE = 0;
ACCEPT_ROOTED_TREES = 1;

Tree tree1 = (anc,der);
tree1.anc.t := 0;

NRM_2 = 0;
MULTIPLY_BY_FREQS_2 = PopulateModelMatrix_2 ("NRM_2");

dynamicFreqs_2 = {4,1};
vectorOfFrequencies_2 = {4,1};
vectorOfFrequencies_2 [0] := computeEQF_2 (AC_2,AT_2,CG_2,CT_2,GT_2,CA_2,GA_2,TA_2,GC_2,TC_2,TG_2);
vectorOfFrequencies_2 [1] := dynamicFreqs_2[1];
vectorOfFrequencies_2 [2] := dynamicFreqs_2[2];
vectorOfFrequencies_2 [3] := dynamicFreqs_2[3];

Model GRMModel_2 = (NRM_2, vectorOfFrequencies_2, MULTIPLY_BY_FREQS_2);
FREQUENCY_SENSITIVE = 0;
ACCEPT_ROOTED_TREES = 1;

Tree tree2 = (anc,der);
tree2.anc.t_2 := 0;

/*---------- Unconstrained Model ----------*/

LikelihoodFunction theLnLik = (filter1, tree1, filter2, tree2);

Optimize (paramValues, theLnLik);
unconstrLik = paramValues[1][0];
unConstrainedParamCount = paramValues[1][1];
fprintf  (stdout, "unconstrained param count: ", unConstrainedParamCount,"\n");

fprintf (stdout, "NRM:\n", NRM, "\n");

/*----------- Constrained Model -----------*/

LikelihoodFunction theLnLik = (filter1, tree1, filter2, tree2);

/* Constraint */
t_2 := t*dynamicFreqs_2[2]/dynamicFreqs[2];

Optimize (paramValues, theLnLik);
constrainedLnLik = paramValues[1][0];
constrainedParamCount = paramValues[1][1];
fprintf (stdout, "constrained param count: ", constrainedParamCount, "\n");

fprintf (stdout, "NRM_2:\n", NRM_2, "\n");

/* Compute and print the log ratio */
LRT = 2 (unconstrLik - constrainedLnLik);
fprintf  (stdout, "LRT: ", LRT, "\n");


/*--------------- function to compute equilibrium frequencies ---------------*/

ffunction computeEQF (dummy1,dummy2,dummy3,dummy4,dummy5,dummy6,dummy7,dummy8,dummy9,dummy10,dummy11)
{
     t = 1;
     tempMatrix = NRM;
     for (k=0; k<4; k=k+1)
     {
           tempMatrix[k][3] = 1;
     }
     tempMatrix = Inverse(tempMatrix);
     for (k=0; k<4; k=k+1)
     {
           dynamicFreqs[k] = tempMatrix[3][k];
     }
     return dynamicFreqs[0];
}

ffunction computeEQF_2 (dummy1,dummy2,dummy3,dummy4,dummy5,dummy6,dummy7,dummy8,dummy9,dummy10,dummy11)
{
     t_2 = 1;
     tempMatrix = NRM_2;
     for (k=0; k<4; k=k+1)
     {
           tempMatrix[k][3] = 1;
     }
     tempMatrix = Inverse(tempMatrix);
     for (k=0; k<4; k=k+1)
     {
           dynamicFreqs_2[k] = tempMatrix[3][k];
     }
     return dynamicFreqs_2[0];
}
[/code]

The two data sets being loaded are modeled data with equal ACGT content (seq1) and double AT content (seq2), where the instantaneous substitution rate to any other base is 1%.

Marissa