Hi,
I am working on a project that involves determines the transition/ transverion ratio of various structures.  We are wondering whether it is possible to get confidence intervals for the estimates of the substitution rates using Hyphy.  
Thanks,
Jen Knies

Dear Jennifer,

I wrote a little example on how to use HyPhy to estimate the ratio and compute the confidence interval using profile likelihood and parametric bootstrap (the latter is more appropriate for smaller datasets). The file can be found at Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login. Download this file and run it in HyPhy to see what it does.

Also, you may find it useful to take a look at pages 9-12 of Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login for some background on confidence intervals and also a walk-through on how to estimate them via the interface.

Cheers,
Sergei

Thanks for the reply.  I downloaded and ran your example file (HKY85_CI.bf) using the p51.nex dataset and tree.  However, I am getting an error message:

Error:
Operation Transpose is not defined for 0.0761426
Current BL Command:Build Formula:meanBS=Transpose(save_bootstrap_kappas["1"])*sa
ve_bootstrap_kappas

This error relates to one of the last lines in the code - computing the mean BS value- and the value (.0761426) is the first value in the save_bootstrap_kappas array.  Can you suggest a way to eliminate this problem?  Thanks!

Dear Jennifer,


wrote on Mar 23^rd, 2006 at 2:52pm:

Error: Operation Transpose is not defined for 0.0761426 Current BL Command:Build Formula:meanBS=Transpose(save_bootstrap_kappas["1"])*sa ve_bootstrap_kappas

This is probably because your version of HyPhy is more than 4-6 months old. The code above uses new matrix indexing options (they are not necessary, but allow for neater looking code), which were added a few months back. Please download the most recent version of HyPhy and try again. Let me know if the error persists.

Cheers,
Sergei

Dear Sergei,
        Thanks- You were right, downloading the latest installment of HYPHY fixed this problem.

I have a question about the output from the HYK85_CI.bf.  I have modified this batch file so that the substitution rate matrix used is the one chosen by data monkey for my dataset (011010).
011010_Rate_Matrix = [*, B*t,t,B*t}
                {B*t,*,B*t,t}
                {t,B*t,*,B*t}      
                        {B*t,t,B*t,*];

In which B replaces kappa_inv everywhere in the batch file.  

The problem is the output values for B here are the inverse of my estimates of B determined by the sample.bf in which the rate matrix has been replaced by 011010.  

Is it correct to assume that the output from HKY_CI.bf is the inverse of the paramater being estimated? And so I should take the inverse of the estimates of "B" (kappa_inv) to have the true value of B.  
Thanks, Jen
???

Dear Jennifer,


wrote on Apr 4^th, 2006 at 6:38am:

I have a question about the output from the HYK85_CI.bf.  I have modified this batch file so that the substitution rate matrix used is the one chosen by data monkey for my dataset (011010). In which B replaces kappa_inv everywhere in the batch file.   The problem is the output values for B here are the inverse of my estimates of B determined by the sample.bf in which the rate matrix has been replaced by 011010.   Is it correct to assume that the output from HKY_CI.bf is the inverse of the paramater being estimated? And so I should take the inverse of the estimates of "B" (kappa_inv) to have the true value of B.

I wrote HKY_CI.bf so that the model is parameterized by the inverse of kappa (i.e. the transversion/transition ratio in place of the more common transition/transversion ratio). This is done for technical (speed of convergence) reasons. However the file 'undoes' the reciprocation when it prints out the output. In your case, you probably don't want to inverse 'B' when it's output. You will notice the last two lines in HKY_CI.bf do the inversion:


Code (]
fprintf (stdout, "\n\nMLE Transition/Transversion estimate: ", 1/kappa_inv, "\nProfile likelihood 95% CI :", Format (1/cmx[2):

, 10, 5), " - ", Format(1/cmx[0], 10,5), "\n"); fprintf (stdout, "\n\nBS  Transition/Transversion mean: ", bs_iterates/meanBS[0], "\nBootstrap 95% CI :", Format (1/upperBS, 10, 5), " - ", Format(1/lowerBS, 10,5), "\n");

.


In your case you should probably write:


Code (]
fprintf (stdout, "\n\nMLE B estimate: ", B, "\nProfile likelihood 95% CI :", Format (cmx[0):

, 10, 5), " - ", Format(cmx[2], 10,5), "\n"); fprintf (stdout, "\n\nBS  B mean: ", meanBS[0]/bs_iterates, "\nBootstrap 95% CI :", Format (lowerBS, 10, 5), " - ", Format(upperBS, 10,5), "\n");

HTH,
Sergei