A simple method for data partitioning based on relative evolutionary rates

RatePartitions

Although it is technically incorrect to use the word ‘partition’ when referring to a data subset, we use ‘partition’ in that sense since this is commonly done in phylogenetics. When partitioning is carried out using TIGER, one must take into account the general properties of the data. One of these properties is that with standard DNA sequence data of protein-coding genes, one to two thirds of the data typically consist of invariable characters. These tend to be binned together to the exclusion of other data when using the TIGER binning strategy or the k-means algorithm (Baca et al., 2017). A partition made of only such data can lead to very high likelihood estimates for it although the assumption that all invariant sites observed in the dataset are truly invariant is most likely wrong. Thus, it is advisable to include a number of slowly evolving characters to create a data partition with low variation. To deal with that problem we developed RatePartitions–an algorithm which works in the following way. The dataset is first run in TIGER or TIGER-rates to calculate the relative rate of evolution for each site (character). These values can range from 1 (invariable sites) to 0 (no common patterns, i.e., the fastest-evolving sites). The sites are then combined into partitions using RatePartitions, which applies a simple formula that ensures a distribution of sites into partitions following the distribution of rates of the characters from the full dataset. This leads to larger partitions for characters with slower rates and, conversely, smaller partitions for those with higher rates. Preliminary tests using MrModeltest v2.3 (Nylander, 2004) and PartitionFinder v.1.0.0 (Lanfear et al., 2012) suggested that this strategy led to models with uniform rate variation within partitions.

RatePartitions is a PYTHON script (Script S1, available at GitHub: https://github.com/jadrankarota/RatePartitions) that determines the rate-spans for a variable number of partitions based on a user-specified division factor and the range of rates calculated by TIGER-rates, and subsequently defines character sets for each partition. The rate-spans are calculated for the first (and slowest) partition with the following function: ${\displaystyle z=x-\left(\left(x-y\right)∕d\right)}$ and for the remaining partitions: ${\displaystyle z=x-\left(\left(x-y\right)∕\left(d+p\ast 0.3\right)\right)}$ where z is the lower limit of the rate-span, x is the upper limit of the rate-span (determined iteratively for each partition, i.e., z becomes x in the following iteration), y is the minimum value of rates for the entire dataset, d is a user defined division factor (which must be greater than 1; a higher number gives a greater number of partitions) and p is the partition number (when >1), which is multiplied by a fixed value of 0.3. The latter reduces the rate-span exponentially as partition number grows, which we found leads to partitions with more uniform rate variation for model-based analyses. Thus, for a dataset with rates ranging from 1 to 0.2 and with d set to 1.5, the first partition will consist of all characters with rates between 1 and 1−((1–0.2)/1.5) = 0.4667. For partition 2, x= 0.4667 and this partition will include characters with rates between 0.4667 and 0.4667−((0.4667–0.2)/(1.5 + 2*0.3) = 0.3397, and so on until less than 10% of all characters are remaining. At this point the iterations are stopped and the remaining characters are placed into their own partition (which becomes the last and fastest-evolving partition).

Simulations

To test the behaviour of our partitioning approach, we simulated DNA sequences on simple trees, then we re-estimated the trees using either TIGER-rates and RatePartitions or the standard partitioning approaches, and we compared the results. We simulated DNA sequences on two different topologies: one asymmetrical (outgroup, P, (O, (N, (M, (L, (K, (J, (I, (H, (G, (F, (E, (D, (C, (B, A))))))))))))))) and one symmetrical (outgroup, (((A, B), (C, D)), ((E, F), (G, H))), (((I, J), (K, L)), ((M, N), (O, P)))). Each tree had 16 terminal branches and all terminal branches were set to a length of 0.1. Using INDELibleV1.03 (Fletcher & Yang, 2009), we simulated for each tree category 14 different datasets of 4,000 base pairs (bp) divided into four partitions of 1,000 bp each. The datasets varied in base composition, substitution models, amount of among-site rate variation, and internal branch length of the tree. We simulated 10 replicates with each set of simulation conditions for both the asymmetrical (AS: simulations on the asymmetrical tree) and symmetrical tree (SS: simulations on the symmetrical tree), resulting in a total of 280 simulated datasets. Full details with all the variables are presented in Table 1. Furthermore, to test for the effect of missing data, we simulated datasets with one of the four partitions randomly removed for three random taxa from the AS1, AS2, SS1, and SS2 sets of simulations (each with 10 replicates). All simulated datasets are available in a public online repository Zenodo (doi: 10.5281/zenodo.1251684), as well as Data S1 in the online version of this article.

From each of the resulting simulated datasets we re-estimated the tree using maximum likelihood implemented in RAxML 8.4 (Stamatakis, 2014) either by (1) partitioning the dataset into the four partitions (4-part) that were simulated, or (2) using TIGER-rates and RatePartitions (TIG). Finally, we compared the resulting trees from these analyses to the original tree used for simulating the dataset using normalized Robinson-Foulds tree distance (RFD) (Robinson & Foulds, 1981). The results are reported as averages across 10 replicates for each set of simulation conditions.

Phylogenetic empirical datasets

We analysed eight previously published lepidopteran datasets (Kodandaramaiah et al., 2010; Sihvonen et al., 2011; Penz, Devries & Wahlberg, 2012; Rota & Wahlberg, 2012; Zahiri et al., 2013; Matos-Maravi et al., 2014; Wahlberg et al., 2014; Rönkä et al., 2016) (Table 2). From the published datasets we excluded sites from the alignment that had more than 80% of missing data unless the whole dataset had 1% or fewer of such sites (data were excluded from Arctiina, Geometridae, Morpho, and Pieridae; final amount of missing data is given in Table S1). All datasets are provided as Data S1 (and deposited in Zenodo doi: 10.5281/zenodo.1252828). The datasets varied in base pair length from 4,435 to 6,372 and in number of taxa from 31 to 164 (Table 2). All datasets included one mitochondrial gene (COI) and four to seven nuclear genes that are commonly used in lepidopteran phylogenetics (CAD, EF-1 α, GAPDH, IDH, MDH, RpS5, wingless) (Wahlberg & Wheat, 2008). We compared 14 partitioning strategies (Table 3), including user-defined ones such as partitioning by gene and by gene and codon position, and a number of different strategies devised based on the relative evolutionary rates assigned by TIGER-rates and division of sites into partitions using the RatePartitions algorithm. We varied the parameter d in the RatePartitions algorithm between 1.5 and 4.5 in increments of 0.5. For comparison of the partitioning strategies we used the BIC score as calculated by PartitionFinder 1.1 (Lanfear et al., 2012). We did two types of searches with PartitionFinder. The first was a user-defined search for direct evaluation of the partitioning strategy obtained with TIGER-rates and RatePartitions. The second was a greedy search (we use “Gr” to mark which strategies included a greedy search), which searches for partitions with similar parameter estimates and combines them so as to reduce the final number of partitions. For example, for a dataset with eight genes that are a priori partitioned by gene and codon position (24 partitions), a greedy search may result in a total of nine partitions because some of the original partitions were combined into a larger subset of data with similar parameter values. BIC was chosen as a statistical model evaluation metric because it usually selects simpler models than AICc and it has been shown to perform well in model selection for phylogenetic analysis (Abdo et al., 2005; Ripplinger & Sullivan, 2008). We refer to analyses with different values of d as TIG1.5, TIG2.0, etc. The greedy search was not performed on TIG1.5, TIG2.0, TIG2.5, and TIG3.0 partitioning strategies because these were shown to have inferior BIC values in preliminary analyses.

We also conducted RAxML analyses of all the eight empirical datasets to compare the resulting trees between those from partitioning by gene and codon position (CodonGr) and the TIGER partitioning strategy with the best BIC score: TIG4.0Gr for Arctiina, Coenonymphina, Geometridae, and Morpho; and TIG4.5Gr for Calisto, Choreutidae, Noctuidae, and Pieridae. We used RFDs to quantify topological differences in resulting trees and we carried out a visual comparison of the trees to determine where the topological differences are located.

Phylogenomic empirical dataset

To test the performance of our partitioning approach on a phylogenomic dataset, we analysed a dataset that includes ca. 1,500 UCEs, nuclear introns from 15 loci, and three mitochondrial gene regions, totalling over 600 kilobases for 18 taxa of gallopheasants and is a difficult dataset because of a presumed rapid evolutionary radiation that occurred in this group (Meiklejohn et al., 2016). Using TIGER-rates and RatePartitions with d = 10, 15, and 20 and running a greedy search in PartitionFinder, we identified the best partitioning strategy as TIG15. Finally, we analysed the dataset using IQ-TREE v.1.6.0 (Nguyen et al., 2015; Chernomor, Von Haeseler & Minh, 2016) with 1000 ultrafast bootstrap replicates (Hoang et al., 2018) either by partitioning according to TIGER-rates and RatePartitions or following Meiklejohn et al. (2016) partitioning strategy.

Simulations

A comparison of trees resulting from analyses of the simulated datasets divided into four partitions (four-part; replicating simulation conditions) shows that they performed slightly better than those with trees resulting from TIGER-rates and RatePartition analyses (TIG) as measured with RFD: across all the simulations on asymmetrical trees (AS) the RFDs for the four-part analyses have a mean of 2.8 ± 4.6, while TIG analyses have a mean of 3.9 ± 6.2; and across all the simulations on symmetrical trees (SS) for the four-part the mean was 4.3 ± 7.1, while for TIG it was 5.7 ± 8.1 (Table 4, Fig. 1). In four datasets, for each type of trees (AS1, AS3, AS4, and AS13; SS1, SS3, SS4, and SS14) the two partitioning strategies systematically recovered the correct tree topology. Overall, more incorrect trees were inferred for the SS datasets than the AS datasets and this was consistent across the two partitioning strategies. The two partitioning strategies behave identically when the difficulty of the dataset increased. Datasets simulated with extremely short internal branches of 0.001 (S5, S6) or with an extreme amount of among-site rate variation together with relatively short branches (S8, S10, S11, S12), resulted in almost all cases in incorrect trees for both the 4-part and TIG analyses. The most extreme departures from the original trees occurred for simulations on trees with extremely short internal branches of 0.001. Under such conditions, even partitioning the dataset according to the initial simulation conditions failed to recover the correct topology, barely performing better than TIGER partitioning (RFD 4-part partitioning: AS5 mean 10.2 ± 3.8, AS6 mean 14.8 ± 6.3, SS5 mean 15.0 ± 4.2, SS6 mean 24.6 ± 3.5; TIGER partitioning: AS5 mean 12.8 ± 4.3, AS6 mean 21.2 ± 6.9, SS5 mean 19.6 ± 5.2, SS6 mean 26.6 ± 1.9). (Table 4, Fig. 1).

Our partitioning method performed equally well even when 25% of the data were missing for three taxa in the alignment (Table 5, Fig. 2). In the AS1 set of simulations, both 4-part and TIG strategies recovered the true tree in all cases, whereas the TIG partitioning performed slightly worse than the four-part in the AS2, SS1, and SS2 simulation sets, as was expected, but RFDs were quite low (RFD 4-part partitioning: AS1 mean 0 ± 0, AS2 mean 0.8 ± 1.0, SS1 mean 0 ± 0, SS2 mean 0.6 ± 1.3; TIGER partitioning: AS1 mean 0 ± 0, AS2 mean 1.0 ± 1.4, SS1 mean 0.4 ± 0.8, SS2 mean 1.0 ± 1.1) (Table 5, Fig. 2).

Phylogenetic empirical datasets

The eight empirical datasets analysed covered a range of taxonomic ranks within Lepidoptera, from genus level (Morpho and Calisto), subtribes (Arctiina and Coenonymphina), two small to medium-sized families (Choreutidae and Pieridae, with about 400 and 1,100 species, respectively), to two very large families (Geometridae and Noctuidae, with over 23,000 and 11,000 species, respectively) (van Nieukerken et al., 2011). The amount of missing data was quite variable. The most complete dataset, Coenonymphina, had more than 90% of sites with less than 20% of missing data, while the least complete dataset, Arctiina, had only 21% of sites with less than 20% of missing data (Table S1).

TIGER partitioning resulted in a different number of partitions for each dataset, with Geometridae and Pieridae being split into many more partitions than the other datasets (Table S2). For example, at d equalling 4.5, Morpho, the dataset with fewest taxa, was split into only seven partitions, Pieridae into 20, Geometridae into 24, while all the other datasets ranged 10–14 in their number of partitions.

In all cases partitioning by gene region was clearly the worst way to subdivide the data, as determined by BIC scores, and applying the greedy search made little improvement (Fig. 3, Table S3). In all datasets, the model with partitioning using TIGER and RatePartitions was the best fit for the data. However, in two datasets (Geometridae and Pieridae), partitioning by gene and codon position with a greedy search came close to the best TIGER strategy, although the BIC scores were still significantly better for the TIGER strategy (Table S3). In all datasets, the improvement in the BIC score from TIG1.5 to TIG3.0 was quite steep, but further differences between TIG3.5, TIG4.0, and TIG4.5, with and without greedy search were relatively small, although the analyses with the greedy search always received a significantly better BIC score. TIG4.5Gr was the best strategy in Calisto, Choreutidae, Noctuidae, and Pieridae, whereas TIG4.0Gr was the best strategy in Arctiina, Coenonymphina, Geometridae, and Morpho (Fig. 3, Table S3).

A comparison of the trees resulting from the RAxML analyses of the best TIGER partitioning strategy with those from the CodonGr analyses revealed that the different datasets responded differently to the two partitioning strategies. For two datasets, Choreutidae and Pieridae, RFDs showed very low values, ca. 5% of internal branches were different between the two partitioning strategies (Table S4). However, most datasets showed great variation between the partitioning strategies, with >20% of internal branches differing. The differing branches tended to be very poorly supported regardless of partitioning strategy (Fig. S1).

An examination of the plots of the relative evolutionary rates estimated by TIGER-rates for each gene fragment and codon position reveals differences among gene fragments, as well as sites belonging to the same codon position in the same gene fragment (Fig. 4, S2). As expected, in general, first and second codon positions receive a much higher rate (i.e., implying slower change) than third codon positions, but in some genes there is a large proportion of third codon positions that also receive a rate of one, e.g., in the Morpho dataset for CAD, EF-1 α, and RpS5 (Fig. S2). Conversely, there are genes that tend to have some fast-changing first and second codon positions, which then receive a relatively low rate. This is usually the case in COI, the mitochondrial gene, but also in several nuclear genes (wingless in all datasets, but also CAD, MDH, and RpS5 in some of the datasets; Fig. 4, S2).

Phylogenomic empirical dataset

When we applied TIGER-rates and RatePartitions on the phylogenomic dataset from Meiklejohn et al. (2016) we recovered an almost identical tree except for the following three differences (Fig. S3). The relationships between P. cristatus, G. gallus and A. rufa shifted from the well supported (P. cristatus, (G. gallus, A. rufa)) in Meiklejohn et al. (2016) to poorly supported (A. rufa, (G. gallus, P. cristatus)). Two other nodes in Meiklejohn et al. (2016) had low bootstrap values (63 and 55). In our analyses using IQTree and the Meiklejohn et al. (2016) partitioning strategy, the bootstrap values for these same two nodes were much higher (93 and 82 respectively), and with the TIGER partitioning the bootstrap for the first node even further increased to 97, while for the second node it decreased to 74. Finally, we observed a great decrease of branch length at the base of the tree when using TIGER.