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.

[StefanFlaumberg]

Hi @benjamin-lieser,

Thanks for the explanations! I will look forward to your preprint to learn more details on how you've managed to overcome overfitting and sensitivity to the guide tree topology. Sounds like a real breakthrough.

Since you're having the preprint almost ready, I guess you would like the code to be merged as soon as possible. Then we will have to fix conflicts with PMSFR and refine the interface/input handling already after your merge. Hope it won't take long.

And yet some more against the -ft AUTO feature:
"is really just a convenience, for inexperienced people wanting a tree in one command."
Sure, the intention is perfectly understandable, but I fear that behind the superficial appeal of simplicity there lie more harm and troubles than utility.
"Our model is very insensitive to the guide tree."
That's great. But, as I wrote earlier, the standard PMSF and PMSR approaches are not, which was already clear from the PMSF original paper. So this option can lure more users into the lazy solution of performing the guide tree search under the simple LG+F+G4 model instead of being mindful about which model might be the best for their data and setup. Since I don't know yet how universally applicable the MutSel model is, I maintain this concern.
"MutSel only works with Amino acids, which I think is also true for PMSF, correct?"
Technically, not correct. Currently, PMSF can be run as well on the DNA and MORPH data types, though uncommonly. You can see in the Issues section someone trying PMSF profile inference from their custom mixture model on MORPH data.
"The parent process waits until the child does the guide tree, so only 2 threads are used with -nt 2"
Yes, only 2 active processes always. But is it strictly guaranteed that the CPUs that ran the parent processes are freed or repurposed to the child processes, i.e. that they are not just stuck waiting? I'm not well-versed with parallelization, but surfing the net shows that technically it is not guaranteed and problems with SLURM can rarely occur.
Also I'm almost sure this won't work properly for iqtree-MPI: at the very least, we need the new process being spawned only by master and setting a barrier for other workers. But sorry if I'm wrong here.

All of these surely are just minor points that could be solved later, after merging the main changes. But I think they are worth mentioning anyway. It's also noteworthy that process spawning, though perfectly fine in general, looks kind of hacky in this context and possibly could be substituted with something like what is done in runMultipleTreeReconstruction(), but this surely would require some work to implement. As all of my initial arguments still stand, in my view, the -ft AUTO option is just an unnecessary code complication with questionable practical value.

Best,
Stefan

[bqminh]

A few comments from Minh and Thomas, which we added also into the source code of the PR. We'll merge it now. But some comments for future release:

Major:

• It doesn't seem to us that you normalise the site-specific Q matrices across all sites. We understand that each individual Q per site is not normalised. But then you still need to make sure that all mu (normalising factor) across sites must normalise to 1, in order to have branch lengths interpreted as the number of substitutions per site. More specifically: (1/nsites) * sum_i mu_i should equal to 1.0, where mu_i = -sum_j pi^i_j * Q^i_{jj}.

Others:

• If possible (for later), rename site_rate_matrices to ptn_rate_matrices
• This is hard coded to contain 190 entries, which only works for AA data and REV model. Better change it to nstates*(nstates-1)/2 for REV and nstates*nstates for NON-REV models.
• What is Cargo.lock file in musel_rust? Is it really necessary to include this file within IQ-TREE code?
• What's the reason for changing cmake req from 10.12 to 11.0.1? Make sure that IQ-TREE can still run on old Mac, such as with Intel chips
• the design of spawning process is not great because you need to load the alignment again causing overhead. It's better to run inside IQ-TREE. Suggest that we do not announce -ft AUTO option for now, and just ask users to explicitly infer a guide tree themselves.
• Better have some option like --max-mutsel-brlen with the default of 200.0. So users can change. We don't like hard-coded constant in the code. Also this has to be in the utils.cpp file, not here, which is hard to debug.

[benjamin-lieser]

This is already a hardcoded constant in IQTREE (10.0), do you want to leave this?

[benjamin-lieser]

The overhead is minimal, but I agree that doing this internal would avoid the arguably a bit ugly process spwan.

I am also open to remove this feature.

[bqminh]

It's not only about spawning a new process, but this actually doesn't work with several other options in IQ-TREE, that you may have overlooked. @StefanFlaumberg pointed out some. Other features like -g option (constraint tree) is not considered (i.e. your guideline won't obey the constraint). A better neat design is to call function runTreeReconstruction(...) instead. You can set params.model_name appropriately, before calling this function.

[benjamin-lieser]

Maybe a general discussion about how to store these values. For the site frequencies it is a vector of pointers, which is arguably brittle with respect to memory leaks or even use after free.
I think a simple vector of vectors should be fine, I doubt we loose performance here.

Would you agree to do this for both?

[bqminh]

I understand. But better leave site frequencies untouched for the moment, as Stefan might touch this variable in another PR. So I prefer not doing both at the moment.

[benjamin-lieser]

A few comments from Minh and Thomas, which we added also into the source code of the PR. We'll merge it now. But some comments for future release:

Major:

* It doesn't seem to us that you normalise the site-specific Q matrices across all sites. We understand that each individual Q per site is not normalised. But then you still need to make sure that all mu (normalising factor) across sites must normalise to 1, in order to have branch lengths interpreted as the number of substitutions per site. More specifically: (1/nsites) * sum_i mu_i should equal to 1.0, where mu_i = -sum_j pi^i_j * Q^i_{jj}.

They are normalized in this way in the Rust part.

Others:

* If possible (for later), rename site_rate_matrices to ptn_rate_matrices

* This is hard coded to contain 190 entries, which only works for AA data and REV model. Better change it to nstates*(nstates-1)/2 for REV and nstates*nstates for NON-REV models.

I agree, I will rework this. But maybe we merge the other PR which touch this first.

* What is Cargo.lock file in musel_rust? Is it really necessary to include this file within IQ-TREE code?

It makes the build 100% reproducible, because all dependencies are pinned. I think this is useful and also how IQ-TREE handles most of it's dependencies (except Eigen and Boost)

* What's the reason for changing cmake req from 10.12 to 11.0.1? Make sure that IQ-TREE can still run on old Mac, such as with Intel chips

10.12 does not support C++17 and is not supported by apple anymore at all. You can run 11.0 on an Intel Mac and one absolutely should. Any hardware from 2013 should support it.

* the design of spawning process is not great because you need to load the alignment again causing overhead. It's better to run inside IQ-TREE. Suggest that we do not announce -ft AUTO option for now, and just ask users to explicitly infer a guide tree themselves.

I have put a comment to the code here.

* Better have some option like --max-mutsel-brlen with the default of 200.0. So users can change. We don't like hard-coded constant in the code. Also this has to be in the utils.cpp file, not here, which is hard to debug.

I also but a comment to the code for this

[bqminh]

I see thanks

[StefanFlaumberg]

Hi @bqminh and @benjamin-lieser,

Since these are directly related to my previous bug-fixing PRs and the PMSR model PR, I address your first two comments in #211. Rate matrices are stored now in ptn_rate_mat similarly to ptn_state_freq, both of them are properly copied by the alignment-copying functions. The added num_rates member of Alignment makes ptn_rate_mat compatible with all data types and non-reversible models for potential usage beyond MutSel. However, I kept the hard-coded 190 constants in the mutsel_wrapper.cpp file, since for now MutSel can only model reversible evolution of protein data.

EDIT:
A comment on the ptn_rate_mat member type (initially, I entirely overlooked your discussion on this topic):
Using vector<double*> for both ptn_rate_mat and ptn_state_freq is not brittle, but rather smart in the context of how things are done in the Alignment::regroupSitePattern() function, while memory leaks or double deletes are ruled out by the design of the Alignment::~Alignment() destructor. Obviously, it is better for ptn_rate_mat and ptn_state_freq to have similar names, types and layouts.

Best,
Stefan