Mutsel implementation

[benjamin-lieser]

This is a implementation of a novel Phylogenetic model, based on Mutation Selection.

It does depend on PhyloGrad, therefore this introduces Rust in the building process.
I also updates the windows and Mac Ci build to use C++17.

It generalizes the PMSF optimization, and allows arbitrary rate matrices. The estimation of the MUTSEL model parameters happens in Rust code, using an automatic differentiation framework.

The final binary does not have any Rust specific dependencies, which allows to distribute static binaries for the supported platforms.

[StefanFlaumberg]

Hi @benjamin-lieser,

A very interesting development!
As I've been recently fixing bugs in the PMSF model and working on the new PMSR model, I've also had thoughts about the possibility of implementing a site-specific matrix model. However, I've ended up concluding that such a model would not be useful in general, as one can reliably neither pre-estimate site-specific matrices in the posterior mean fashion nor use site-specific matrices with free parameters (such as most DNA models have). Both approaches suffer from the same problem of overfitting.
But since your approach invokes PhyloGrad, which somehow seems to circumvent this overfitting problem, and produces site-specific matrices with fixed parameters, this seems to be exactly the case where site-specific matrix model usage is justified.

Surely, there will be some code conflicts between the changes in this PR and in my PR for the PMSFR model. Most notably, in the way the ModelSet object is constructed. There are also some problems here regarding how site-specific matrices are handled by different functions of the Alignment class (now they are not handled at all, yet they should be). And, possibly, the file-reading functions could be merged or be using a common helper.
Someone will have to resolve the conflicts anyway. In my plans, I was going to make the PMSFR model PR fully ready for merge (finally) in the coming two weeks. If you agree to have your PR be merged after mine, then, after the PMSFR model PR is accepted, I could send a PR with the conflict-resolving changes directly to your branch so that you could check that everything is alright before the changes are sent here.
Otherwise, I'll just have to fix the conflicts alongside updating the PMSFR model for the current release (that code is one year old), which is also acceptable, but might take more effort to be done correctly.
What feels more convenient for you?

The only thing that seems very off to me in your changes is the -ft AUTO option spawning a new process to run iqtree in it for a full tree reconstruction. There are a number of serious drawbacks to this approach:

1. While the option looks very general, using the LG+F+G4 model makes it compatible only with the protein data. This is a latent bug.
2. Using site-homogeneous models to infer guide trees is known to result in the PMSF profiles being more susceptible to LBA artifacts. It is quite possible that using a frequency mixture model, such as LG+C20+G4, even with the default mixture weights and --fast or -n 10 can produce more reliable PMSF profiles in roughly the same time. Actually, this is quite an important research direction!
3. The option itself (even if given a better implementation) is unnecessary in practice. Every computational biologist is totally capable of manually performing the 2 consecutive iqtree runs required by the PMSF approach. As discussed above, it is better to have to perform the first, guide-tree-inferring run without shortcuts, explicitly controlling the accuracy-speed tradeoff by choosing the right model and the search stop rule.
4. Adding a new header file only to implement an optional feature is a clear overkill. The iqtree codebase is already very large and full of dead code, so we should be rational about such additions.
5. The new-process-spawning approach is totally new to the iqtree codebase and for a good reason: having an iqtree run calling another run seems cumbersome given how complicated the program is. Apart from this abstract architectural argument, my concerns are about multi-threading and MPI. For instance, if I run with -nt 2 -ft AUTO, am I guaranteed to be using only 2 CPUs? This is quite important for running on HPC clusters. I cannot emphasize more that the whole thing with all its complexities and concerns is completely unnecessary here.

Take a look at this from the following angle: the new code, possibly not parallelization-safe, is added only to allow for a feature that is not only unnecessary, but also promotes a potentially suboptimal practice!
Therefore, I strongly recommend this feature be removed.

Best regards,
Stefan

[benjamin-lieser]

Hi,

thanks for the comments. We want to have a preprint out soon, so there the details of the model will be explained. Maybe two key points here. Overfitting is not prevented by Phylograd, but by regularizing the parameters against a posterior mean model (PMSF like, but made a bit faster). The second is that Mutsel allows us to use no rate categories at all, all the rate variation can be explained extremely well just by Mutation Selection (as is would make a lot of sense biologically).

Having a common file format to specify site specific models would probably be very useful. I guess the features needed would be exchangebility, equlibirum and rate for every column (or I guess pattern, we do not need different parameters for the same pattern). Plus an optional RHAS model. So if you implement such a thing in the PMSFR model, we could merge after and use this.

The -ft AUTO option is really just a convenience, for inexperienced people wanting a tree in one command. Our model is very insensitive to the guide tree.
I can see that the self spawning logic is quite a bit (mostly because it is very robust), but I am pretty sure most of the people using IQTREE just want simply the best practice algorithm for their tree and a guide tree could push them to use Model Finder over PMSF/Mutsel which will for sure yield worse results.
MutSel only works with Amino acids, which I think is also true for PMSF, correct?
The parent process waits until the child does the guide tree, so only 2 threads are used with -nt 2, in any case you do not need to use it and I guess it is very unlikely someone has the guidetree file called AUTO.