Ultrafast classical phylogenetic method beats large protein language models on variant effect prediction

Thank you for your thoughtful review. Please find our response below:

The idea of simplifying the trees by partitioning the MSA into pairs of similar sequences is an approximate representation of the cherries in the tree. Since the partitioning algorithm is a greedy one, it may overlook the full structure of the MSA, reducing the relationships between sequences to binary relations between pairs. That also explains why CherryML with FastCherry lose the accuracy compared to CherryML with FastTree in Figure 1b.

The whole purpose of developing FastCherries is to avoid doing expensive computations while assuring little reduction in accuracy. It is a worthwhile tradeoff that can benefit applications requiring scalability, such as the one considered in our manuscript. Indeed, for the variant effect prediction task, it can be seen in Supplementary Figure S1 that FastTree is too slow to scale up to the full dataset size. In contrast, FastCherries provides no loss in performance and can crunch all the available data, yielding the best results.

This work represents an incremental improvement over the original CherryML proposal, and therefore, its novelty is somewhat limited.

We respectfully disagree with this view. As our results clearly demonstrate (Figure 1, Supplementary Figure S1, and Supplementary Table S1) and as the other reviewers have noted, the performance of CherryML with FastCherries is comparable to that of CherryML with FastTree while being one to two orders of magnitude faster. Furthermore, we are proposing a novel method (SiteRM) to estimate site-specific rate matrices and this framework has the potential to advance phylogenetic inference significantly; software such as IQTree and its partition model puts this application well within reach. When compared to the seminal LG model, SiteRM provides a substantial improvement in variant effect prediction, as seen in Supplementary Table S1, suggesting that it can model site-specific evolutionary constraints more accurately than previous approaches.

For the machine learning community, it is not sufficient to demonstrate that one method is better than another through benchmark results alone. A deeper analysis explaining why the phylogenetic model performs better than language models in certain benchmarks and when it fails to do so would provide valuable insights for the community to learn from.

How could the authors provide a deeper analysis of why the phylogenetic model performs better than language models in certain benchmarks?

The improved performance of our method comes from conditioning on the wildtype to obtain a local fitness function $P(y | x, t=1).$ Please see our response to Reviewer 2 for more discussion on the intuitions. We plan to make this more clear in the final version. We do agree with the reviewer that deeper error analysis of large language models would be a valuable research direction. Ultimately, we believe that the best variant effect predictors will combine our evolutionary modeling $P(y | x, t)$ with protein language modeling.

Could you please explain a bit on why on Figure 1.b when the number of family increases the error rate seem increase for CherryML?

The fact that the (rounded) error rate is 1.7 for $2^{11}$ and 1.8 for $2^{13}$ is just a small sampling effect. At these larger sample sizes, the statistical risk of the estimator has converged to its bias term.

Thank you for your thoughtful review. We respond below:

The paper lacks a comparison to other methods in terms of computational resources, such as memory usage. Could the author provide some results?

Our method uses linear space; the logarithmic factors in the computational runtime come from divide-and-conquer and binary searches. This makes our method space efficient (up to a constant). We will discuss this and add practical results on memory usage in the final version. In more detail, the space complexity is $O(s^2 br + nl + ls).$ Here, $O(s^2 br)$ is the cost of the precomputed matrix exponentials, which is shared across all MSA’s. $O(nl)$ is the cost to load one MSA into memory. With a runtime of $n \log n$, pairing takes $O(n)$ space, since it takes $O(n)$ to store the distance between x and all other sequences, and the space complexity equals the cost to store the stack of recursive calls. $O(ls + r)$ is the space of computing the initial site rates. The $O(ls)$ term comes from the count matrix $\text{count}[i,j] =$ the number of times character $j$ occurs at position $i$. Computing site rates given branch lengths takes $O(r)$ space and computing branch lengths given site rates takes $O(l)$ space. In both cases, we only track the best rate for the current site or best branch length for the current site.

The model is still highly parameterized. Can the SiteRM model avoid overfitting with extensive datasets?

As described in Section 3.3 of our manuscript, we avoid overfitting by using pseudocounts drawn from the LG model, which serves as a prior. For the results described in the paper, we find that mixing the pseudocounts with the data at a 1:1 ratio works well in practice, but it should be possible to tune the regularization coefficient $\lambda$ using cross-validation. Other options we did not explore that would benefit from large datasets would be some form of weight-sharing of the rate matrices across different positions.

How does the method perform with datasets that have a large number of gaps in MSAs (for example, the sequence has very low homology)?

The QMaker "insect" dataset is especially notorious for containing a large number of gaps. We observe a similar performance to FastTree in Supplementary Figure S1.

Can the methods be extended to incorporate other types of biological data, such as DNA/RNA sequence or structural alignment?

Indeed, by using compound states one can model the joint evolution of amino acids and structural states. Our software is not limited to amino acid alphabets and can be applied to learn site-specific rate matrices for other kinds of molecular data such as DNA/RNA or structural states. We plan to pursue this in future work.

Thank you for your thoughtful review. Please see our response below:

The authors state "The first probabilistic model of protein evolution was proposed by Whelan and Goldman". It would strengthen the paper to better cite the models that Whelan and Goldman were inspired by.

Thank you. We will cite the earlier models in the final version, in particular Dayhoff's and JTT's.

I would appreciate the author comparing to more state-of-the-art mutation effect predictors.

In Supplementary Table S1, we show results including all variant effect predictors from the ProteinGym paper.

How does this work relate to the conclusions of Weinstein, Eli N., Alan N. Amin, Jonathan Frazer, and Debora S. Marks. 2022. “Non-Identifiability and the Blessings of Misspecification in Models of Molecular Fitness and Phylogeny.” Advances in Neural Information Processing Systems, December.?

While we agree with the theoretical results presented in the above paper, our work challenges the assumption in their theorems that there exists a unique, "global" fitness function f. Instead, our intuition is that the fitness function is contextual, i.e. it is different in different parts of the evolutionary tree. For example, in some subtree, a positively charged amino acid may be required at position $i$, while in another subtree a negatively charged one may be required. This may be because of different contexts the protein is evolving in, due to e.g. protein-protein interactions. Our approach to variant effect prediction captures this intuition by using the "local" fitness function $P(y | x, t=1)$ at each point in the tree. We do not think that the theoretical results in the referenced paper necessarily explain our strong performance on variant effect prediction. Instead, we believe that it is the use of "local" fitness functions that is leading to our improved performance. There is certainly exciting theoretical and empirical work to be done in this direction.

How do you parameterize the substitution matrix?

We use the parameterization from the CherryML paper, namely the off-diagonal elements of $Q$ are given by $\text{diag}(1/\sqrt{\pi}) \cdot S \cdot \text{diag}(\sqrt{\pi})$ where $\pi = \text{softmax}(\theta)$ for an unconstrained vector $\theta$ (here $\pi$ will be the stationary distribution of $Q$), and $S$ is a symmetric matrix given by $S = \text{SoftPlus}(\Theta + \Theta^T)$, where $\Theta$ is an unconstrained upper triangular matrix. This unconstrained parameterization allows CherryML to use unconstrained first-order optimizers to quickly estimate $Q$, as implemented by libraries such as PyTorch and Tensorflow.

Thank you for your detailed review. We respond below:

while the paper demonstrates the effectiveness of SiteRM in variant effect prediction, further exploration of its applicability to other evolutionary biology tasks or datasets could further understand its capabilities and limitations.

We agree and plan to explore other applications of SiteRM such as improving phylogenetic tree inference in future work. Software such as IQTree and its partition model puts this application within reach.

there lacks more detailed information on the benchmarking process, including datasets used, and potential biases in the evaluation, which could improve the transparency and reproducibility of the results.

We focussed on standard benchmarks from prior work (CherryML and ProteinGym) with the hope that it maximizes transparency and reproducibility. We also provided code with detailed instructions to reproduce all our results. Nonetheless, we agree that our work is not self-contained and more details on the benchmarks could be included. We plan to incorporate these more detailed descriptions of the benchmarks into the final version.

it would be better to provide a more vivid and more understandable method to explain the provided algorithm, such as drawing some diagrams or flowcharts of the pipeline.

This is a great suggestion. The closest we currently have to this is the runtime analysis in Appendix A.4.2. We will include a more graphical depiction of the algorithm in the final version.

although the provided method can improve the end-to-end runtime of the process, its performance drops a little comparing with the original method. Could you please list some possible reasons to explain the performance decrease, and some possible solutions.

FastCherries shows a small asymptotic bias due to the use of Hamming Distance (HD) in the pairing step. We explored other alternatives during the project which we decided not to include in the paper, such as (1) pairing based on maximizing the composite likelihood of the pairing, (2) pairing based on minimizing the MLE distance between pairs, and (3) even a random pairer. The random pairer actually worked quite well for the WAG model of protein evolution on some benchmarks, showing asymptotic consistency with a relative statistical efficiency of ~1/8th, but worked poorly for the LG model due to the challenges of estimating site rates. We found that approaches (1) and (2) can indeed provide more accurate estimates in some cases. Since pairing based on HD worked very well already, we decided to make it the focus of the paper. We plan to explore other pairing methods in future work. It would be exciting to find a variant of FastCherries that retains the near-linear runtime while showing an even smaller error.

the performance of the method may vary across different protein families or evolutionary contexts. It would be better to assess the generalizability of the approach to diverse datasets and evolutionary scenarios.

Our benchmark on the QMaker datasets (Supplementary Figure S1) shows that our method performs well on datasets from diverse parts of life. We agree that more in-depth error analysis may provide further insights and improvements to the method, but we leave this for future research.

the accuracy of the estimated rate matrices and predictions might depend on the quality of the input multiple sequence alignments. Are there any ways to address potential biases or errors in the MSAs to improve the robustness?

This is an important question which pertains to much work done in statistical phylogenetics. Although we did not explore it in our work, the speed of our method makes it possible (unlike prior work) to obtain bootstrap confidence intervals for the rate matrix estimation, which should enable users to understand the extent to which the estimates are stable to e.g. subsampling the MSA (whether rows or columns) or changing the MSA building algorithm.

the benchmark used to evaluate the performance may have limitations or biases. It would be better to compare the performance on other benchmarks.

As mentioned above, the QMaker benchmark provides evidence in this direction showing generalization to diverse domains of life.