Hi Sergei,

I have been trying to compute the likelihood for a codon model and have found that the answer I get with HyPhy is different from the one I get if I compute the likelihood "by hand"; that is, by multiplying together the transition probbilities and summing over all ancestoral nodes.

I have attached the bf file for a simple (fake) data set (seq.fasta) comprising of 3 taxa with only a single observed codon site which is the same for all taxa. The rate matrices and equilibrium frequencies, which were previously estimated and are contained in the file ParameterEstimates.txt, are allowed to differ between branches as indicated in the specification of the tree topology.

Please could you help me understand why my computed likelihood (as given in ComputeLikelihood.bf) differs from that computed by HyPhy.

Much appreciated!

Miguel

Thanks so much for your quick reply! Your code works with this example, but I can't seem to get it to work with a more complicated example (probably due to an unrelated problem!)...

I have attached another batch file (ComputingLikelihood3.bf) which is essentially an extension of the previous one I sent you. Here I have expressed the EFVs r0, r1 and rx as functions of two variables k0 and k1 as well as a previously estimated EFV (pi.txt). An estimated rate matrix is given (Q.txt) but must be multiplied by the equilibrium frequencies, so you end up with three different rate matrices and EFVs which apply to the different branches of the topology as specified in the Tree statement. The simplified topology consists of 3 taxa with only one observed codon site for each taxon. Taxon N1 has codon AAA and taxa N2 and N3 each have codon GGA.

Given Q, pi, k0 and k1, I have computed the likelihood of the tree. But, if I compute the likelihood a second time by simply "copy-and-paste"-ing the previous code below it, I get a different likelihood value! This is relevant as I intend looping over (k0, k1) pairs and computing the likelihood for each pair, which currently does not work.

Also, neither of the above two likelihoods equal the likelihood that I have computed manually, which is also concerning.

Your help would be greatly appreciated!!

Thanks in advance,
Miguel

Dear Sergei,

I am trying to implement a GY94 model of codon substitution and cannot understand one of the results I'm getting.

I have attached a simplified version of the problem where most of the parameters (alpha, kappa, mu and t) are fixed. The only parameter (k) which is estimated through optimization of the LF occurs in the EFV.

The frequency r_{i} of codon i is computed as:

[code]
global k = 0.5;
k:>-0.000000000000000000000001; k:<1.00000000000000000000000000001;

r={Rows(pi),1};
        
for(i=0; i<Rows(pi); i=i+1){
   
   if(pi[i][1]==wt){
     /* codon codes for wild AA */
     global r[i][0]:=k*pi__[i__][0]/aa_freq__[wt__][0];
   }
   else{
     /* codon does NOT code for wild AA */
     global r[i][0]:=(1-k)*pi__[i__][0]/(1-aa_freq__[wt__][0]);            
   }
}
[/code]

where pi is the "usual" EFV computed using HarvestFrequencies, aa_freq is a vector containing the frequencies of the amino acids and wt is an index corresponding to the "wild type" amino acid occurring at a given site.

The code I have given you focuses on only a single site in Gag where the only amino acid that is observed across all HIV strains is E (so that wt=16 corresponding to E). In such a case, I would expect the MLE of k to be 1 (or at least very close to it). However, the MLE according to HyPhy is very close to 0 - the exact opposite of what I would expect! This would mean that the columns in Q corresponding to the codons encoding E are all zero vectors - surely this cannot be the case!?

I have played around a lot with this and have discovered that if I make all the branch lengths equal to one, the MLE for k is 1 as expected - I'm not sure why though.

Hope you can help. Thanks a lot!

Miguel

This is really weird. I was using the development version you told me to previously in this thread and get:

Code:

Likelihood Function's Current Value = -34.5934741151561
k=0.00383036
mu13.3622==13.3622
kappa9.78031==9.78031
alpha0.653415==0.653415 

I have now downloaded the latest version and get the same answer! What could I be doing wrong?

Miguel

Hi Sergei,

I'm sorry for being stupid  , but where can I get hold of the version of HyPhy that you are using to get the correct answer above? I have downloaded the latest hyphy-r298.zip from Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login, but I'm still getting the wrong answer.

Thanks again,

Miguel

Thanks! This makes sense!!

Perhaps you could also help me understand the following observation: Whenever I set all the branch lengths to a very small number, say 0.01, k goes to zero, but setting all the branch lengths to a larger number, say 1, produces the desired result k=1. I can't quite grasp why this should be the case.

Cheers,
Miguel