Reviewer 2: Robert Beiko
(The review by Robert Beiko includes two consecutive rounds of comments, here interleaved. The second exchange is indented.)
In any event, the paper is certainly worthy of publication — I have some continuing problems with some details that I outline below, but this is a great contribution to an ongoing discussion.
This paper is a detailed response to Koonin and Wolf (ref [2
]) in which a compositional model was used to generate phylogeny-free sequences to test the model selection framework of ref [1
]. The paper consists of (1) A justification of the need for a formal test of common ancestry, (2) A detailed description of several critical limitations of null hypothesis-based tests of statistical significance, and (3) a detailed refutation of the work of Koonin and Wolf.
I think part (1) is fine; I have not followed the genetic code debate closely, but I agree that quantitative hypothesis tests based on evidence apart from the genetic code are worthwhile.
Part (2) addresses some well-covered areas of statistical theory, although the "Wolf vs. Koonin" etc. example is less illuminating than the rest of the section. I think e-values get short shrift, though: if one considers two different e-values from two different comparisons, I think one would be justified in taking the larger of the two e-values to constitute greater evidence for a lack of process generating similarity between the two sequences (convergence, common ancestry, etc.).
Response: This is an exceedingly common misinterpretation of frequentist null-hypothesis significance tests. As explained in the text, by the very logic of null-hypothesis tests, a large E-value (or P-value) is evidence for nothing. This cautionary point is routinely reiterated in introductory statistical texts, e.g. "Although a 'significant' departure provides some degree of evidence against a null hypothesis, it is important to realize that a 'nonsignificant' departure does not provide positive evidence in favour of that hypothesis. The situation is rather that we have failed to find strong evidence against the null hypothesis."
]. Similarly, "… a large significance level indicates that, with respect to the particular test used, the data provide no evidence that this hypothesis is false. This should not be interpreted as evidence in support of the hypothesis, but merely as a lack of evidence against it. "
Here is one way of understanding why it is fallacious to make conclusions from insignificant results. If the null hypothesis is true, then we expect to find probable results from the null distribution (i.e., we expect a large "insignificant" E-value, e.g. E > 0.1). But using probable results under the null to conclude that the null is true is committing the fallacy of affirming the consequent: if A then B, we observe B, therefore A — an invalid syllogism. In the case of sequence similarity tests, observing low similarity (e.g., a small S statistic in BLAST search) is expected under the null distribution of random alignments. But we also expect low similarity if the sequences are homologous and have diverged for a large amount of time.
Given the lack of well-developed alternative models to homology, it seems reasonable to me to consider a small BLAST e-value as constituting stronger evidence for homology than a larger one. I agree that this is not as strong as explicit model comparisons, but it does nonetheless highlight an important problem with the entire procedure as reported in the 2010 Nature paper and here.
Before I get to part (3), I have some major concerns about the approach, in particular how the tests are formulated and how the results are interpreted as positive evidence of common ancestry.
First and foremost, any likelihood or Bayesian model-contrasting framework is only as good as the models you feed to it.
Response: This is true, but from a Bayesian point of view, this is not really a problem. Rather, it is a fundamental and inherent property of the scientific method itself, common to all methods of inference in general, and Bayesian methods handle it in the only proper way. We never know the true model, and in fact we can be certain that even our best models are wrong. To quote Box: "all models are wrong, but some models are useful"
]. The purpose of a Bayesian analysis (and of science in general, in my opinion) is this: given a set of models hypothesized to explain our data, which model in this set should we believe is best at explaining the data? The most we can ever do is to say that one model is more probable than another (and quantify how much more probable)
The original 2010 paper contains the sentence "Statistically significant sequence similarity can arise from factors other than common ancestry, such as convergent evolution due to selection, structural constraints on sequence identity, mutation bias, chance, or artefact manufacture". But none of these alternatives is actually tested alongside the hypothesis of common ancestry!
Response: I did indeed test some of those alternatives. In my original 2010 Nature paper I explicitly tested models of independent ancestry and chance convergence of the observed sequence similarity (as well as the nested, hierarchical patterns in the real protein data). Furthermore, my independent origin models are good approximations of independent origin models that include convergence due to several of those other factors (see below).
The independent origin tests use amino acid substitution models that were generated under the assumption of homology, and apply them to sequences that are completely homologous or very nearly so (see next paragraph). We seem to all be in agreement that these sequences are homologous (indeed the author says so in his response to Yonezawa and Hasegawa [116,120]: "neither selection nor physical constraints alone can plausibly generate the high levels of sequence similarity (55% average sequence identity) observed in the universal protein data set that I used."
Response: My quoted claim is not equivalent to conceding that the sequences are homologous. There are other explanations, besides homology, for high sequence similarity. My statement was not an assumption but a conclusion based on my model selection results, taking into account realistic biophysical and biochemical considerations.
This statement taken in isolation is fine, but is at odds with the claim that a meaningful, quantitative hypothesis test has been formulated and carried out, since no attempt has been made to develop plausible alternative models (for instance, models with independent origins + convergence) and let them fight it out with competing alternatives.
Response: First, I did formulate and carry out a test with alternative models —
Figures and from the original paper
]are graphical representations of some of these alternative, independent ancestry models. The independent ancestry models that I considered could in principle produce the observed protein data, but the model selection tests show that the probability of chance convergence from independent origins is extremely low
Second, my response to Yonezawa and Hasegawa
]provides an explanation of why my independent ancestry models are reasonable approximations to models explicitly incorporating convergence due to physical constraints and selection for function. Any independent origin model incorporating selection and/or physical constraints must be based on the known biophysical and biochemical properties of proteins. The most pertinent characteristic of proteins is their extreme functional redundancy. To a good first approximation, primary sequence determines structure, and in turn structure determines function. Yet many sequences give the same structure, and many structures provide the same function
]. All our biophysical and biochemical data indicate that the sequence space consistent with any particular function is so large as to be nearly indistinguishable from random background residue frequencies
]. For selection to converge on similar sequences, sequence space must be tightly coupled to function space, so that independent selection pathways for the same function would be likely to result in similar sequences. In reality, however, there is an incredibly large sequence space associated with any particular function — so large that knowledge of the function alone tells us essentially nothing about the possible sequences that could specify that function
In accord with these biophysical facts, each independent ancestry in my models has a prior on the residue frequencies given by their background frequencies (these can be provided a priori by the residue substitution matrices, but in my analyses they were actually inferred from the protein sequence data). Hence my independent ancestry models approximately account for (1) constraints on protein sequences imposed by function and (2) the prior probability of selection for function converging on similar sequences. For the same reasons, the independent ancestry models can also be considered as representing design hypotheses in which a "designer" intended to produce proteins with similar functions and/or structures, with an ignorance prior on the sequence space fulfilling those criteria (reflecting the fact that we know nothing about the intentions of the designer). For proteins with highly dissimilar sequences, it would be important to carefully and explicitly account for the very weak constraints on sequence imposed by function. But the proteins in my universal data set have relatively high sequence similarities. Hence the extremely large model selection scores strongly suggest that these approximations are reasonable, since subtle differences in background frequencies will have a negligible effect on the model selection scores reported in my Nature letter.
The "universal" proteins do indeed have high sequence similarities, but at the same time the ribosomal proteins you examined would seem to be an extreme case in comparison with the enzymes that were examined in e.g., [126
]. The highly interacting nature of ribosomal proteins, along with their tendency to co-localize in large operons, would potentially make the probability of convergence quite a bit higher than for your typical garden-variety "active sites + a bunch of glue" enzyme. I agree that a proper model of independent origin of (at least some) ribosomal proteins + convergence would not be expected to tip the scales in terms of model evaluation, but the degree of model preference could easily shrink by quite a bit.
Response: This is a good point — physical constraints, such as necessary physical interactions between ribosomal proteins and the ribosome, increase the probability of convergence (i.e., increasing constraints decrease the size of the relevant sequence and structure space consistent with a given function). For proteins with very weak similarities, these constraints may be important to model. In any case, my results are robust to omission of the ribosomal proteins.
I mention "completely homologous" because there is a serious problem contained within at least some of the supplementary examples of the 2010 paper. Consider for instance the three sequences in Section 4.3, where E+M have significant sequence similarity, as do M+P, but not E+P. But looking at the BLAST alignments clearly shows these sequences for what they are: M is a multidomain protein that is partially homologous with E and partially homologous with P in a way that does not overlap with E. But the alignment columns are then shuffled, such that M+P are identical in columns 1,2,5,8,9, and so on, while E+M are identical at 3,4,6,7, etc. Were these sites properly clustered, inspection of the BLAST local alignments would immediately give a clear indication of what is going on. But since the columns are shuffled (and I know of no process that would give rise to this pattern of intercalated similarities), BLAST is railroaded into showing both local alignments as covering the entire protein length. This seems like an irrelevant example when the highly conserved proteins used in the main paper will show no such tendencies. I would like to know how the example in 4.2 was generated as well.
Response: None of this matters for the intended point — these sequences were not intended to mimic proteins in my dataset, they were made to illustrate a methodological point. The point of the example in Section 4.3 is simply that the model selection tests prefer independent ancestry for certain highly similar sets of sequences (sequences with extremely statistically significant similarity, E 10-100). Therefore the methodology considers factors other than mere sequence similarity in assessing independent vs common ancestry, contra Koonin and Wolf
The second problem is that the entire notion of common ancestry is ill-defined, especially in light of lateral genetic transfer. I think it is unlikely in the extreme that the components of the universal cellular machinery examined here are anything but homologous, but it seems to me that establishing common ancestry of each of these components does not guarantee that there was a single ancestral organism (or population of organisms, for any reasonable definition of 'population') in which all copies of the machinery were ancestral to the ones we would see today, which would seem to be necessary to make the leap from "common ancestry of a bunch of proteins" to "a common ancestor of life". This matter has been raised by Peter Gogarten, Olga Zhaxybayeva, Ford Doolittle, and Joel Velasco (pers. commun.), among others. What exactly is being formally tested here?
Response: Universal common ancestry is defined here as "the proposition that all extant life is genetically related"
]. The fundamental question I addressed was not whether the universally conserved genes in the Brown dataset are similar — they unquestionably are. Rather, I tested whether these sequences — and by necessity the organisms that carry them in their genomes — are genetically related by common ancestry. If you carry genes in your genome that are homologous to genes carried in the genomes of all other living organisms, then logically you are genetically related to all other living organisms. It is not simply that we have "common ancestry of a bunch of proteins". Rather, we have "common ancestry of a bunch of fundamental proteins that are each found in all known living organisms"
Lateral genetic transfer is an orthogonal question. All life could be genetically related, via only vertical genetic inheritance from ancestor to descendant. Alternatively, all life could be genetically related, with rampant lateral gene transfer. There are many intermediate possibilities.
The issue of the nature of the "universal common ancestor" is also largely independent of the status of universal common ancestry. For instance, Ford Doolittle speaks of common ancestry without a common ancestor. It will be useful to quote him at length:
We (some of us) do doubt that there ever was a single universal common ancestor (a last universal common ancestor or LUCA), if by that is meant a single cell whose genome harboured predecessors of all the genes to be found in all the genomes of all cells alive today. But this does not mean that life lacks 'universal common ancestry' — no more than the fact that mitochondrial DNA and Y-chromosome phylogenies do not trace back to a single conjugal couple named Eve and Adam whose loins bore all the genes we humans share today means that members of Homo sapiens lack common ancestry. That 'common ancestry' does not entail a 'common ancestor' is perhaps a subtle point …
Similarly, I wrote in the original 2010 Nature letter:
UCA does not demand that the last universal common ancestor was a single organism, in accord with the traditional evolutionary view that common ancestors of species are groups, not individuals. Rather, the last universal common ancestor may have comprised a population of organisms with different genotypes that lived in different places at different times.
The fact of the matter is that organisms existed in the past that contained ancestral copies of these fundamental proteins, and that all known extant life has inherited descendant versions of these proteins from these ancestral organisms (at least that is the scenario that my model selection tests strongly support). Whether you want to call that ancestral group of organisms "a species", "a population", "an ancestor", "the LUCA", or something else is primarily an issue of semantics. My universal common ancestry models are also consistent with a single individual organism/cell that contained all the ancestral versions of these proteins, though that is extremely unlikely given basic pop-gen considerations (and even less likely with LGT).
Because all extant living organisms carry descendants of these fundamental proteins, at some point the ancestral group of organisms necessarily had the ability to exchange fundamental genetic material. This is one argument for calling this ancestral population a "species", regardless of whether it was exchanging genetic information via sex or LGT or whatever — e.g., that would be a strict application of the Biological Species Concept. Similarly, Doolittle and others have argued that "the proper way to model prokaryotic evolution over 4 Gyr is as a single, albeit highly structured, recombining population, not an asexual clade."
The rest of this review I dedicate to the Results and Discussion, part (3) of the paper. 3a, 3b and 3c refer to the three claims that are to be refuted (top of p. 15), 3a is further split into four parts.
3a-I, pp. 15-17: Of course phylogeny of a group of sequences models the hierarchical structure of a set of sequences in a way that simple pairwise BLAST cannot. But sequence similarity is evidently a major contributor to the optimal tree in terms of its branch lengths and topology, and based on the examples presented in the SI of the original Nature paper, it seems that one could use a graph-based approach to building putatively homologous sequence sets from cliques or almost-cliques of BLAST results. To me, it seems that extending the analysis of sequence similarity to sets (n > 2) of sequences may be the critical factor in the model test, rather than the following step of phylogenetic inference and modelling. Can a counterexample be constructed where most or all sequences in a set show statistically significant similarity to one another, yet the likelihood-based tests support independent ancestry over common ancestry?
Response: Yes, and I presented exactly this type of counter-example in Section 4.3 of the Supplementary Information of my original paper. I have also constructed many other examples where three or four sequences all have highly significant sequence similarity to each other (rather than just two of three pairwise comparisons as in Section 4.3). In response I now present an additional four-sequence example in the text.
3a-II, p. 17: Similarly to above, I don't think it's fair to say that the three proteins in the example in 4.3 of the Nature paper show "significant similarity" when only 66% of the sequences show such similarity.
Response: Homology is transitive — if A and B are homologous, and B and C are homologous, then A is homologous to C. If significant sequence similarity is in fact a reliable marker for homology, then those three proteins should be homologous even though only two of three pairwise comparisons show significant sequence similarity. Analogous examples are commonly found in real proteins. Regardless, the new four-sequence example should mollify this concern.
Homology is clearly not transitive in the case of multidomain proteins. For example, if protein P1 has domains X+Y and protein P2 has domains Y+Z, it is not correct to say that a third protein P3 containing only domain Z is homologous to P1 because P2 happens to contain elements of both P1 and P3.
And this is clearly what is happening in the three- and four-protein examples: careful inspection of the proteins in Appendix 2, Additional File 1
shows that each column either supports a perfect partitioning of the sequence sets into BE+MP, or BM+EP, or BP+ME. So this example is equivalent to taking six distinct, non-homologous domains, and building sequences E = D1D3D5, M = D1D4D6, B = D2D3D6, P = D2D4D5 with odd/even pairs forced to sit on top of one another in an alignment. Shuffling the columns to generate the Appendix 2, Additional File 1
sequences allows BLAST to construct local alignments that span non-homologous "gaps" between the fragments of homologous domains, thus making it seem like the sequences are similar to each other in non-exclusive ways.
My point here is that these examples have similarity patterns that are evolutionarily implausible,
Response: This is exactly the point. The patterns of similarity in these artificial sequences are highly unlikely to have been generated by evolutionary, phylogenetic processes as described by the common ancestry models. These patterns of similarity, however, could easily result if the sequences have independent origins (e.g., from some sort of convergent process or by design — the latter independent origin hypothesis is the correct one for these sequences). Whatever the case, the sequences are highly similar, and yet nevertheless the model selection tests (correctly) choose independent ancestry over common ancestry.
and the argument that model comparisons capture more than simple sequence similarity are not bolstered by these artificial scenarios. The manner in which these examples were constructed could be described as "partial common ancestry", and raise the question of what model should be preferred when sequences are forcibly aligned in such a way that part of the alignment is true homology and part is not. The preference for BE+MP in Table is presumably a consequence of more columns supporting this pairing than any other.
The 4.2 example is more confusing to me as I do not see how the likelihood is generated for the independent origin model. Is it based solely on stationary amino acid frequencies, without any branch lengths estimated? It seems to me that an example with two sequences is not ideal to illustrate the distinction between the two types of approach: were a third sequence added that was similar to 1&2 such that all pairs of the triad were significantly similar, I presume that common ancestry would once again be supported.
Response: I agree that the 4.2 example is less applicable, so I have removed reference to it in the text.
3a-III and IV, pp. 18-19: I do not think K&W suggested that the universal proteins could have been generated by a process lacking in phylogenetic structure:
Response: As I read it, K&W implied that their profile model was a suitable model for the universal proteins, at least as far as my model selection tests were concerned. Otherwise, what is the relevance of their profile model to the actual data analyzed in my paper?
rather, the point of their paper was that the key result of CA > IA could be generated from a dataset in which CA was not the correct model.
Response: As I argued in detail in the text, K&W have never shown that their artificial dataset was not generated by a CA model.
It also seems odd to contrast a star tree model with the partial independent origins model, when the natural complement to a star tree would be an independent origin for each of the sequences considered. I expect under those circumstances that the star tree would be favoured over a true independent origins model, as is the case with the K&W data.
Response: I also expect that a star-tree would be preferred over independent origins for each of the sequences, but I don't find that to be a very informative comparison. In contrast, the star-tree vs independent ancestry tells us something important. Star-trees model similarity, but not nested hierarchical patterns. The independent origin models do the opposite; they model hierarchical patterns well, but not similarity between the independent groups. The key observation is that star trees are worse than independent origin models for the real proteins (as gauged by the model selection scores). This strongly suggests, then, that hierarchical structure is more important than similarity for the real protein data.
What I was getting at here is that in the partial independent origins model, hierarchical structure and similarity are still conflated at every level except that of the deepest interdomain relationships. And I would expect these relationships to have the smallest overall impact on the question of hierarchy vs. similarity, since they would be the least stressed out about having to coexist in a star tree, as opposed to much more similar sequences that desperately want to share an internal branch somewhere, anywhere. I don't think the distinction here is as clean as it is made out to be.
Response: I agree the distinction is not completely clean. Nevertheless, a hypothesis that can describe some hierarchical structure, but absolutely no similarity between two groups ("partial" independent ancestry), does much better than a hypothesis that cannot model hierarchical structure, yet models similarity between the groups well (star tree). The independent ancestry hypothesis cannot account for any similarity between the two groups, whereas the star tree can, and yet the star tree, a common ancestry model, does much worse. Therefore, something in the independent ancestry hypothesis, something necessarily other than mere similarity, is responsible for the independent ancestry hypothesis beating the star tree. The only factor that could be responsible is hierarchical structure, and this is the case even though only the deepest interdomain relationships are really affected.
Concerning the profile model, it's not surprising that the model under which the data were generated turns out to perform best. Clearly the K&W simulation was not intended to mimic the structure of the universal protein data set, but the fact remains that data sampled under residue probability distributions, with no true phylogenetic structure, generated a data set which led to an incorrect preference for the model presented in the original paper.
Response: As I explained in the text, if the data are generated by a common ancestry model (as K&W's profile model is), and the tests choose common ancestry, then they chose correctly. And more importantly, when the true generating model is included in the set of models being tested, the model selection tests choose the true generating model.
Also, please clarify why the amino acid dataset has a data-to-parameter ratio of 0.63.
Response: For the profile model, each column in the alignment has 19 free parameters (the 20 residue frequencies that must be estimated, which sum to one). In the real, universal protein data set, there are only 12 sequences, so each column has only 12 observations from which to estimate the 19 parameters, giving a data-to-parameter ratio of 12/19 = 0.63.
3b: Under "K&W's justification is faulty", two definitions of "independent" are conflated. What I take from the K&W quotation is that the patterns of similarity at each site are independent from one another, since there is no underlying tree. Conversely, the common assumption of "independently evolving sites" does not require that sites have uncorrelated histories or patterns of similarity, rather that changes at one site are not conditional on changes at other sites.
Response: Admittedly the K&W quotation is somewhat vague as to just what they mean by "independent". However, it does not matter for my point if your interpretation is correct. Just because the patterns of similarity at each site are independent from one another does not mean that there is "no signal of common ancestry" or "no evolutionary signal", as K&W claim. As I mentioned in the text, a star tree produces patterns of similarity at each site that are independent from one another — and a star tree is a bona fide phylogeny, being an evolutionary representation of common ancestry.
K&W themselves never claimed that their sequences were artificially convergent, only that they are highly similar.
Response: Three counter-examples from their comment: "these alignments contain no signal of common ancestry (in more general terms, no evolutionary signal) whatsoever … ". "Alignments of statistically similar but phylogenetically unrelated sequences successfully mimic the purported effect of common origin." "The computational experiment on the estimation the likelihoods of alignments of similar but evolutionarily unrelated sequences." Sequences that are similar yet evolutionarily (or phylogenetically) unrelated are by definition convergent at the level of sequence
While these sequences were indeed sampled from homologous alignment columns, they certainly would possess what I would consider to be "conflicting phylogenetic signal". Indeed there seems to be some confusion about whether that term is meant to refer to reticulate evolutionary patterns within a single set of homologous genes, or the partial homology problem I point out above.
Response: By "phylogenetic signal", I mean non-random, nested hierarchical patterns of similarities in the sequences. It is possible for a single set of sequences to have two or more very different, incompatible nested patterns of similarities, which is what I mean by "conflicting phylogenetic signal". This can occur in concatenated alignments of homologous proteins when there is reticulate evolution (since each protein is treated as one single long sequence). But it can also occur within a single protein (for instance reticulate evolution of domains in a multi-domain protein), or within a single protein domain. For my purposes the cause of conflicting phylogenetic signal is irrelevant — the point is that conflicting phylogenetic signal is highly unlikely to result from a common ancestry model (in which a local region of sequence is modelled with a phylogenetic tree). Regardless, any apparent conflicting phylogenetic signal found in K&W's artificial data is necessarily only due to random chance.
The comparison of the sampling protocol of K&W to the models of Bruno is very interesting, although the fact that the generating columns had hierarchical structure to begin with makes me wonder whether the analogy is accurate or not. In any case, does this mean that we cannot in fact distinguish homologous proteins from nonhomologous proteins that exist in a reduced but shared sequence space due to structural constraints?
Response: In some cases, I think the answer is yes. However, imagine we have a reduced sequence space and yet the actual sequences show strong nested, hierarchical patterns. In that case the common ancestry models should explain the data better than independent ancestry models (better than models that account for the reduced sequence space yet cannot account for hierarchical patterns), and so we should be able to distinguish the two.
3c: You should probably add the assumption that the sequences used to generate the substitution models (JTT, WAG, etc.) were themselves homologous, and emphasize the assumption that the Markovian process is stationary, reversible, and homogeneous.
Response: I've modified the text regarding the Markovian assumptions. But the substitution models deserve some comment, and I have added a section of clarification in the text. First, for the formal purposes of the model selection tests, we don't actually know that the substitution models (like WAG) were based on homologous sequences, and this assumption is unnecessary. The only thing we really know is that the sequences used in constructing WAG, etc. were similar. From a model-selection perspective, the origin of the substitution matrices is irrelevant — the likelihood function doesn't know where the matrices came from, and the likelihood is the final judge. All that matters is that some substitution matrices give higher likelihoods than others.
But the WAG matrix was built on the underlying assumption that a phylogenetic tree related the sequences whose changes contributed to the entries in the matrix. If WAG was in fact forcing similar but non-homologous sequences to obey a phylogenetic framework, wouldn't that have some impact on your contrasting of models that have varying components of phylogenetic relatedness?
Response: The WAG matrix is, ultimately, nothing more than a proposed set of substitution rates. Regardless of where those rates came from, they may work well in the evolutionary models or they may not (relative to other rate matrices). More importantly, perhaps, is the fact that we can dispense with an assumed rate matrix altogether and use Mr Bayes to infer a general time reversible amino acid rate matrix from the universal protein data itself. I have done such with the universal proteins; it doesn't change the results. And the matrix inferred from the actual data is qualitatively quite close to WAG.
These assumptions seem particularly relevant if we are testing common ancestry of sequences that may have initially emerged in a period of reduced amino acid diversity and other drastic differences in evolutionary processes. I'm not an expert on early amino acid alphabets, but it seems that contrasting CA vs. IA at a time when there were 10 or fewer amino acids in use would drastically change the degree to which different models are preferred, especially if a meaningful model for convergence could be developed.
Response: These considerations are worth keeping in mind. But since all of these proteins (from widely divergent taxa) have the same 20 amino acids and the same genetic code, the most parsimonious explanation is that the genetic code (and amino acid alphabet) was set in its modern form well before the divergence of this set of proteins (however, Syvanen offers an alternative view that arginine and tryptophan were added rather late