Customizing the constraint on the omega parameter of YN98

[Keren Halabi]

Hi dear team,

I am experiencing an issue with a branch-site codon model I've written in Bio++.

The model consists of 2 mixture model, where each of these consists of 3 sub-models of type YN98 each.

Each mixture model consists of 3 omega categories, such that: omega_0 <= omega_1 <= 1 <= omega_2

each of these omegas is mapped to an omega_ parameter of a sub-model of type YN98.

The omegas in the mixture models are parameterized by the following parameters:

p in (0,1)

omega1 in (0,1)

omega2 in (1,999)

k in (0,50)

such that:

omega_0 = (omega1 * p) ^ k

omega_1 = omega1 ^ k

omega_2 = omega_2 ^ k

The issue:

Each omega_ parameter in the corresponding YN98 model is bound within [0.001,999]. This means that:

(lower_bound(p) * lower_bound(omega1) ) ^ upper_bound(k) must be greater of equal to 0.001

upper_bound(omega2) ^ upper_bound(k) must be less or equal to 999.

These two restrictions are limiting the bounds I can set for p, k and omega2 more than I would like. Since that function that allows to set the constraint on the YN98 model parameters is protected, I am unable to customize the constraint on the omega_ parameters from the constructor of my model.  Thus, I was wondering if there is any way I could add a public function to YN98.cpp that would allow it? 

If yes, are there any additional parameters mapped to YN98's omega I should set the constraints on as well? 

If not, how would you recommend to tackle this issue?

Many thanks in advance!

Keren

[Laurent Guéguen]

Dear Keren,

it is not possible, because not safe. Those bounds have been set to prevent the generator to be singular, and then issues
in the computation of the exponentials.

In your code, to prevent exceptions, it is possible to use AutoParameters,
where values outside constraint range are set to the limit (bpp-core/src/Bpp/Numeric/AutoParameter.cpp: line 77).

More on a modelling point of view, it does not make sense to ask for omega to be 0.000001 instead of 0.001 (and the same with high values),
the accuracy of the model is far below this level of precision. Moreover, the biological interpretation of differences between extreme values of
omega (and of most parameters, actually) is questionable.

Cheers,
Laurent


Btw, I notice on the forum that some of your former requests are without reply, but I cannot remember if we settled them by email.
Are you still waiting for a reply to these?

Simulating under a codon site-model Par Keren Halabi 1 message 2 vues 6 sept.

SubstitutionModelSet - changing the node IDs associated with a model after its creation Par Keren Halabi 1 message 3 vues

[Keren Halabi]

Thank you, Laurent!

This sounds exactly like what I need. Additionally, your input on the reason for using the YN98 bounds to begin with made us rethink our parameterization, so it was a very good input.

 

As for the other referenced posts:

Simulating under a codon site-model- Resolved

SubstitutionModelSet - changing the node IDs associated with a model after its creation - Pending for your response.

Many thanks again!

Keren

[Keren Halabi]

Dear Laurent,

I am experiencing some optimization issues due to the lower bound on omega, which is 0.001 in Bio++. Specifically, the omega parameter that represents purifying selection (omega0) often reaches the value 0.001, which makes me believe that if I had the option to further extend the bound on this parameter, I could have potentially achieved a better fit. In such case, setting my parameter as auto parameter would not help, because the value cannot eventually reach a lower point than 0.001.

I see that in PAML, the lower bound on omega0 is 0.0001 (please see line 1852 in codeml.c, which is attached). Is there any way to do the same in Bio++?

Many thanks!

Keren

[Laurent Guéguen]

Dear Keren,

it is possible if you change the limits in the source code in YN98.cpp.

The 0.001 limit is due to computation problems when omega=0, but 0.0001 would do the same job.

However, I would say that, given the gross approximation due to modeling, the biologic interpretation is the same for

omega=0.001 and omega=0.0001, so I do not see why bothering. 

Cheers,

Laurent