A question regarding partitioned phylogenetic analysis of morphological data

[HS6986]

Dear IQ-TREE Developers,

This method of partitioned phylogenetic analysis of morphological data and the command below don't work well. See the replies below.

I am reanalyzing morphological datasets from published studies using maximum likelihood (ML), Bayesian, and maximum parsimony methods. For the ML analysis, I am using IQ-TREE v2.3.6 and trying to partition data by the number of states of each character, as done in MrBayes (see MrBayes Manual, p. 142) as it may increase model fit. This method seems to have been applied to IQ-TREE, for example, in Černý & Simanoff (2023).

Characters recognized as invariant by IQ-TREE in the dataset were removed using Mesquite and my custom R script in order to be able to apply an ascertainment bias correction. A partition file for the dataset was generated by my custom R script. I ran IQ-TREE with the following command:
iqtree2 -s MatrixVarchars.nex -p MatrixVarchars.nex.partitions.nex -m MFP+MERGE+ASC -B 1000 -bnni -alrt 1000 -T AUTO

The execution appeared to finish without any apparent issues, but my question is, whether IQ-TREE can automatically apply the number of states of each partition to each Mk model. I scrutinized the log file and the iqtree file, but I could not find any related information.

This is a zip file of the folder for my reanalysis of Wang et al. (2023). I would appreciate it if you could answer this question.

[bqminh]

That's a good question. For MORPH data type IQ-TREE doesn't print the number of states in the model. But internally, the states are always encoded from 0, 1, 2,.... 9, A, B, ..., Z. (http://www.iqtree.org/doc/Substitution-Models#binary-and-morphological-models). The number of states will be determined by the largest state (in this order) in a partition. For example, if a partition has states 0, 1, 2; the number of states will be 3. If some states are missing, they will be technically removed from the model. For example, if a partition has states 0, 1, 3, 4, 6 (i.e., 2 and 5 do not occur in the alignment), then it will be considered as a five-state model, including only states that occur at least once in the alignment.

Does that answer your question?

[HS6986]

Thank you!

This is the answer I was expecting. Now I can apply IQ-TREE to morphological data with confidence.

Thank you again for this information.

[HS6986]

I just realized that the command above probably doesn't work well.

If you want to partition a morphological dataset, you have to split it into multiple files by the number of states. If you don't do this, IQ-TREE will probably assigns missing (?), ambiguous (or polymorphic) (e.g., {1,2}), and inapplicable (-) characters in each partition of a morphological dataset equal likelihoods (see the FAQ documentation) corresponding to the number of states in the entire dataset, rather than in each partition, which unexpectedly makes k of the Mkv model in each partition equal to the number of states in the entire dataset in most cases (morphological datasets almost always contain a large number of these characters). This was inferred from information in log files of partitioned phylogenetic analysis of other morphological datasets and the IQ-TREE specifications.

I hope that there will be no unlucky people like me who have to analyze morphological datasets again later.

[davidcerny]

Hi @HS6986, David Černý here. The problem actually goes further than that: even if your matrix was entirely free of missing data, ambiguities, and polymorphisms, IQ-TREE would still determine the state space of every character based on the largest observed number of states. For example, even if your dataset contains a single 7-state character but mostly binary characters otherwise, the evolution of all those binary characters will be modeled using a 7-by-7 Q-matrix, and states 2, 3, 4, 5, and 6 will be treated as present in the model but never observed. Because the Mk model has equal stationary frequencies and exchangeabilities, this is an extremely unlikely outcome, so your likelihoods will take a huge hit: see Mulvey et al. (2024). Also, your branch length estimates may become biased: see Khakurel et al. (2024).

Because IQ-TREE considers binary and multistate characters to be different data types, and because it doesn't allow you to mix data of different types in the same file, the only workaround is – as you correctly state – to split your matrix into multiple files. In Černý & Simonoff (2023), we provided an R script to do that; you can find it here. I wrote a slightly more robust version of that script for a class I taught at UChicago back in 2023; you can find the updated version here. I also wrote a blog post a few months ago where I go into detail on how to implement partitioning by state space size in various phylogenetic software including IQ-TREE, RAxML-NG, and RevBayes:

https://davidcerny.github.io/post/teaching_revbayes/

[HS6986]

Dear Dr. Černý,

Thank you very much for providing this information. In fact, the reason I realized that the above method for partitioned morphological phylogenetic analysis didn't work well internally was partly because I took a close look at the data in your paper (Černý & Simanoff, 2023). Thank you for this great work.

The fact that, even if the partitions are specified in a partition file, if a morphological matrix is not split into multiple files, IQ-TREE determines the state space of every character based on the largest observed number of states in the entire matrix, is new to me. Thank you for providing this information. However, your statement that IQ-TREE doesn't accept mixed data of binary and multistate characters in the same file seems odd to me, because, in IQ-TREE, binary characters in morphological matrices can be treated as two-state MORPH characters (by explicitly specifying -st MORPH) and modeled in two-state Mk models rather than in JC2 models. In an IQ-TREE Issue I opened (#410), I provided several examples of analyses demonstrating that this is the case.

I've written an R script to remove characters recognized as invariant in IQ-TREE in a morphological matrix and split it into multiple files by the number of states, using the R packages TreeTools and phangorn, for my studies not yet published. Because of my lack of knowledge of programming itself, R, and its packages, I'm sure there are many improvements and alternative ways for each step, but I hope this helps. I also provide an R script for generating a matrix without characters recognized as invariant in RAxML-NG and a partition file for RAxML-NG (models can be selected by ModelFinder in IQ-TREE, for example). Here's the folder containing these scripts. If you notice any problems with the implementation of them, please do not hesitate to let me know.

I've just started studying mainly phylogeny of Cambrian metazoans, and I've seen several studies in this filed that performed unpartitioned phylogenetic analyses of morphological datasets in IQ-TREE (though always combined with MrBayes and TNT) (e.g., Howard et al., 2020a, b; Shi et al., 2022; Wang et al., 2023, 2024; in this field, ordered (additive) characters seem to be almost never used). I'm concerned that their morphological data were most likely inadequately modeled.

Thank you again for providing the information.

[davidcerny]

Hmm, that's a good point – for 2-state characters, -st MORPH -m MK and -st BIN -m JC2 should in theory specify the same model, and the former has the advantage of allowing such data to be intermixed with multistate characters. This made me question whether the information I'd given you, i.e.,

even if the partitions are specified in a partition file, if a morphological matrix is not split into multiple files, IQ-TREE determines the state space of every character based on the largest observed number of states in the entire matrix

was always true. To test this, I've conducted an experiment using the dataset of Fonseca et al. (2024), which I'm currently using for an unrelated project. First, I split it by character type and state space size into separate Phylip files called MK2.phy, MK3.phy, MK4.phy, MK5.phy, MK7.phy, ORDERED3.phy, ORDERED4.phy, and ORDERED5.phy. I then used the following partitions.nex file:

#nexus
begin sets;
    charset part1 = MK2.phy;
    charset part2 = MK3.phy;
    charset part3 = MK4.phy;
    charset part4 = MK5.phy;
    charset part5 = MK7.phy;
    charset part6 = ORDERED3.phy;
    charset part7 = ORDERED4.phy;
    charset part8 = ORDERED5.phy;
    charpartition mine =
      JC2+ASC+G: part1,
      MK+ASC+G: part2,
      MK+ASC+G: part3,
      MK+ASC+G: part4,
      MK+ASC: part5,
      ORDERED+ASC+G: part6,
      ORDERED+ASC+G: part7,
      ORDERED+ASC+G: part8;
end;

and conducted an "Analysis 1" as follows:

~/iqtree-2.3.5-macOS/build/iqtree2 -p partitions.nex --runs 100 -T AUTO -allnni -safe --prefix analysis1

Then, I combined these files back together to obtain a single Phylip file called morph.phy, which stores a reorganized version of the original matrix where the characters that share the same number of states form contiguous blocks. I then created a corresponding partition file:

#nexus
begin sets;
    charset binary = 1-767;
    charset state3unord = 768-876;
    charset state4unord = 877-884;
    charset state5unord = 885-887;
    charset state7unord = 888;
    charset state3ord = 889-937;
    charset state4ord = 938-941;
    charset state5ord = 942-943;
    charpartition mine =
      MK+ASC+G: binary,
      MK+ASC+G: state3unord,
      MK+ASC+G: state4unord,
      MK+ASC+G: state5unord,
      MK+ASC: state7unord,
      ORDERED+ASC+G: state3ord,
      ORDERED+ASC+G: state4ord,
      ORDERED+ASC+G: state5ord;
end;

and conducted an "Analysis 2" as follows:

~/iqtree-2.3.5-macOS/build/iqtree2 -st MORPH -s morph.phy -p partitions.nex --runs 100 -T AUTO -allnni -safe --prefix analysis2

As it turns out, Analysis 1 found ML trees with log likelihoods of around −25,900, while Analysis 2 only found trees with log likelihoods closer to −38,500. The huge difference is clearly due to the fact that rate matrices of different dimensions were being applied to the data. Unfortunately, this demonstrates that my original claim was correct: IQ-TREE doesn't determine rate matrix dimensions based on the largest number of states in a given partition, but rather the largest number of states in a given file. The only way to force IQ-TREE to employ rate matrices of the right dimension is to split the dataset into multiple files – partitioning it by other means is not enough.

P.S. I'd be interested in your experience applying ModelFinder to morphological data. To me, it never seemed super useful, except maybe for choosing the rate heterogeneity model. My impression is that in morphology, many aspects of model choice are either determined by the background knowledge of the researcher (e.g., when to use MK vs. ORDERED) or they can be straightforwardly derived from the properties of the dataset (when to use +ASC). You could allow IQ-TREE to tell you that, for example, GTR2 fits better than MK, but the fact that character state labels are inherently arbitrary for morphological data makes it very difficult to assign a meaningful interpretation to the resulting estimates of rates 0<->1, 1<->2, 0<->2, etc.

[HS6986]

Dear Dr. Černý,

Thank you for you careful analysis and information.

The IQ-TREE specification you practically demonstrated—even if the partitions are specified in a partition file, if a morphological matrix is not split into multiple files, it determines the state space of every character based on the largest observed number of states in the entire matrix—is very tricky. Thank you for showing this practically. I am convinced. I wish it worked at least as well as RAxML-NG. If I have time, I will try to encourage the IQ-TREE developers to make IQ-TREE behave similarly to RAxML-NG, or to BEAST 2, in this regard (as you can see from #410, I think IQ-TREE has almost never envisioned this kind of use).

I did state that models of morphological data can be selected using ModelFinder, but this only meant selecting rate heterogeneity across sites (characters) (RHAS) parameters (+G4 and +Rx) on MK+ASC (or ORDERED+ASC) by specifying -m MFP+ASC (see #410 for examples of practical use). I apologize for the confusion. I think that assigning +G4 to each partition of morphological data, as you did, is not too suboptimal to be honest, but I believe that selecting models in ModelFinder improves model fit a bit. However, IQ-TREE does not seem to allow you to specify whether a partition is unordered or ordered in advance, for example, in a partition file, and then select the RHAS parameters. Thus, if a morphological dataset has both unordered and ordered characters, I think you have to run ModelFinder on unordered and ordered partitions separately first, although this is not an elegant way and does not allow you to optimize parameters on the entire dataset.

Also, as you correctly pointed out, using GTR2 (and GTRX, which is actually built internally in IQ-TREE) for morphological data is very problematic. By default, IQ-TREE always considers only MK when selecting models for morphological data. We should always use only MK or ORDERED for morphological data.