Allow different mutation rates to different bases while using per-position mutation rates

[Duncan Ralph]

At the moment we're using the feature where we pass in per-position mutation rates for our sequences, e.g.

position, rate
0, 0.3
1, 0.2
[...]
37, 0.1

But we also need to be able to specify different probabilities of mutating to each base, say position 1 almost always goes to C. We'd like to pass in something like this instead:

position, rate_A, rate_C, rate_G, rate_T
0, 0.1, 0.7, 0.1, 0.1
1, 0.85, 0.05, 0.05, 0.05
....
35, 0.02, 0.02, 0.02, 0.94

So two things: 1) we're pretty sure this isn't present in bppseqgen at the moment, but if it is present hopefully someone can tell us? and 2) we're assuming it's not there, so we're about to implement it, so any pointers would be much appreciated.

thanks!

[Laurent Guéguen]

Hi Duncan,

in the case of bppseqgen, it would mean describe a model per site at most.
Actually it is possible, but on a git branch of bpp-phyl that is still under development, called newlik.

In this branch, the computation of the likelihoods is totally rebuilt, and then it is possible to
perform much more exotic things, mixing trees, data and models in as many ways  as possible.

It is still under big development, but if you are interested to take the risk, you are welcome.

Cheers,
Laurent

[Erick Matsen]

Thank you much for your quick response, Laurent.

Is `newlik` destined to become the primary likelihood computation engine?

[Laurent Guéguen]

Yes it will. I do not know when exactly, but it will.

Just now bppseqgen is not updated for this, but if you need it in the next days I can do it easily.

Cheers
L

[Erick Matsen]

Laurent--

It would be really wonderful if you could do that. We aren't in a huge rush-- within the next month would be just fine.

As far as I can tell, the best we can do to "watch" your development is to clone the repository and then pull regularly. Is that correct? 

Thank you,

Erick

[Julien Yann Dutheil]

Hi there,

Just one more note though: you might also need a substitution model that has as parameters a probability of mutating for each state... I am quite sure such a model has been proposed in the literature, but I don't think we have implemented it in Bio++. It would not be difficult to add it, but this would have to be done as well if I'm correct. One model I can think of is to have a JC69 exchangeability matrix with free equilibrium frequencies (that is, a HKY85 with kappa==1), but I am not a 100% sure that would be equivalent to what you are asking for?

Cheers,

Julien.

[Erick Matsen]

Hello Julien--

Yes, a model with free equilibrium frequencies would work great for us. The model that you describe is I believe the F81 model.

A cursory scan of your documentation doesn't show an F81 implementation, but as you say we could use HKY85 with kappa set to 1.

Thank you!

Erick

[Laurent Guéguen]

Hello,

this is done, in newlik branch (both in bpp-phyl and bppsuite repositories).

In newlik, there is a new syntax for config files, to get with the new features.
You have to define  the trees, models,  processes that mix both and simulation procedures
that put processes along the alignment. 

I have attached an example, with command as usual

bppseqgen param=seq.bpp

Every remark and idea is welcome 
cheers
Laurent

ps : since I am a nice guy, I have added F81 model is bio++.

[Duncan Ralph]

Excellent Laurent, thanks! I'll give it a try and let you know how it goes.

[Erick Matsen]

Laurent--

You are indeed a nice guy! We owe you many beers, bottles of wine, etc.-- feel free to cash in next time I'm in your area or you're in my area.

We're currently prepping for an upcoming talk and getting a paper back to reviewers, but this will be at the top of the list after that. Duncan will definitely be stress-testing this feature!

Erick