DiscML: an R package for estimating evolutionary rates of discrete characters using maximum likelihood

Background The study of discrete characters is crucial for the understanding of evolutionary processes. Even though great advances have been made in the analysis of nucleotide sequences, computer programs for non-DNA discrete characters are often dedicated to specific analyses and lack flexibility. Discrete characters often have different transition rate matrices, variable rates among sites and sometimes contain unobservable states. To obtain the ability to accurately estimate a variety of discrete characters, programs with sophisticated methodologies and flexible settings are desired. Results DiscML performs maximum likelihood estimation for evolutionary rates of discrete characters on a provided phylogeny with the options that correct for unobservable data, rate variations, and unknown prior root probabilities from the empirical data. It gives users options to customize the instantaneous transition rate matrices, or to choose pre-determined matrices from models such as birth-and-death (BD), birth-death-and-innovation (BDI), equal rates (ER), symmetric (SYM), general time-reversible (GTR) and all rates different (ARD). Moreover, we show application examples of DiscML on gene family data and on intron presence/absence data. Conclusion DiscML was developed as a unified R program for estimating evolutionary rates of discrete characters with no restriction on the number of character states, and with flexibility to use different transition models. DiscML is ideal for the analyses of binary (1s/0s) patterns, multi-gene families, and multistate discrete morphological characteristics. Electronic supplementary material The online version of this article (doi:10.1186/1471-2105-15-320) contains supplementary material, which is available to authorized users.

Some of these programs have advanced (or realistic) features that are not implemented in other programs. For instance, the BayesTraits program implements adistribution for rate variation [8]. The GLOOME program allows the estimation of prior root probabilities of the character states [10,15]. The BadiRate program allows variable birth rates and death rates, and corrects for unobservable data [13]. Furthermore, many multistate characters do not necessarily evolve in a BD manner [16], and should therefore be modeled using transition rate matrices other than BD.
In order to perform accurate rate estimation on a variety of discrete characters, we have developed a unified http://www.biomedcentral.com/1471-2105 /15/320 program DiscML by implementing the advanced features mentioned above as well as flexible options for transition rate matrices.

Implementation
DiscML estimates the evolutionary rates of discrete characters by fitting the distribution of all character states (the data) on a given phylogeny. The data need to be in a matrix format (vector format for a single site) as required in many other phylogenetic programs in R (see examples in Additional file 1). The provided phylogeny is required to have branch lengths, as branch lengths will be used as a relative time scale in the analysis. The evolutionary rates, transition rate matrices, and additional parameters discussed below will be optimized to maximize the likelihood of the data. The optimization is achieved using the PORT routines [17] implemented in the nlminb function in R.

Implementation of rate variation in the analysis
Rate variation among the character sites has long been recognized and implemented in DNA analyses [18], but has been missing from most analyses of non-DNA discrete characters (but see [8]). DiscML considers rate variation among the character sites by implementing a discrete distribution (with the option of alpha=TRUE).

Estimation of prior root probabilities
Most programs for the analysis of discrete characters assume only uniformly distributed prior root probabilities, e.g., π 1 = π 2 = .. = π a = 1 a , (a is the total number of character states). DiscML allows the estimation of prior root probabilities (π a ) for different character states (with the option of rootprobability=TRUE).

Flexibility on both the transition model and the number of character states
DiscML is flexible on both the size and type of the transition rate matrix (Q), which can be customized by users. This option could open the door for novel evolutionary analyses on different discrete characters. Several transition rate matrices are pre-determined in DiscML: model="ER" (equal rates, i.e., all entries in equation 1 are equal), model="SYM" (symmetric, i.e., α 1 = α 2 , β 1 = β 2 , γ 1 = γ 2 , ..), and model="ARD" (all rates different, i.e., all entries are free to vary). ER and SYM are reversible matrices, while ARD matrices are irreversible.
Finally, all transition rate matrices (Qs) are calibrated [19], i.e., each Q satisfies so that the evolutionary rate parameter (μ) is the average number of transition events per site per evolutionary time unit [20].

Forced reversibility and flexible irreversible options
When the prior root probabilities (π) for different character states are estimated, reversible transition matrices will no longer necessarily result in reversible evolutionary processes (because of potentially different probabilities of character states). Since it is sometimes of biological interest to assume reversibility (i.e., the expected x → y changes equal to the y → x changes), DiscML can allow forced reversibility by setting reversible=TRUE. In practice, reversibility is obtained by multiplying the corresponding root probabilities (equation 4) to the entries in reversible transition matrices, e.g., ER and SYM. Such a practice is conceptually the same with the general timereversible (GTR) DNA substitution model [21]. In Dis-cML, model="GTR" is equivalent to the combination of model="SYM" and reversible=TRUE. Similarly, when the prior root probabilities for different character states are estimated, forced reversibility can be applied to the BD related matrices (equation 5).
In DiscML, the default setting is reversible=FALSE and users have the flexibility to conduct analysis by assuming irreversible evolutionary processes. Unlike in reversible processes, the root position can greatly affect the maximum likelihood calculation in irreversible cases [22,23]. Therefore, it is only meaningful to perform irreversible analysis on a rooted tree. If the provided phylogenetic tree is unrooted, DiscML will first reroot the tree by midpoint rooting, and perform analysis on the midpoint rooted tree.

Correction for unobservable data
Some characters may contain unobservable character states, which can only be inferred indirectly from the presence of observable states of the same characters in related taxa. Ancient characters can be lost from all examined extant taxa, and result in unobservable data. Dis-cML provides the option of zerocorrection=TRUE to calculate the likelihood conditional on a pattern being observable following [24], i.e., where L − is the likelihood of unobservable patterns. The correction for unobservable data (shown as '+0' in Table 1) is essential for systems such as gene family data due to the complete loss of some ancient genes, but not suitable for single-site analyses and for systems in which all character states are observable (e.g., nucleotide bases).

Site and branch specific estimations
Even though the default setting of DiscML is to perform rate estimation by fitting the distribution pattern of all character sites on a phylogeny, there is an option to perform rate estimation on individual sites (ind=TRUE). Individual rates can be graphically displayed using plotmu=TRUE. Furthermore, DiscML allows branch specific rate estimation, which can be specified using '$' on branches in the provided tree file. For instance, (((taxon1$1: 0.01, taxon2$1: 0.01)$3: 0.01, taxon3$2: 0.02)$3: 0.01, taxon4$2: 0.03); specifies three rates, one for the branches leading to taxon1 and taxon2 ($1), one for the branches leading to taxon3 and taxon4 ($2), and one for The parameter μ is the estimated evolutionary rate of the characters. "1s/0s only" indicates binary analysis by converting all non-zero characters to 1s using simplify=TRUE, '+0' indicates the correction for unobservable data using zerocorrection=TRUE, '+ ' indicates the implementation of a discrete distribution using alpha=TRUE, '+π ' indicates the estimation of prior root probabilities using rootprobability=TRUE, '+π REV ' indicates the estimation of prior root probabilities with forced reversibility using rootprobability=TRUE and reversible=TRUE.
the remaining branches ($3). The modified tree files are no longer in the conventional Newick format, we have developed a function read.tree2 in DiscML to read such modified tree files.

Additional features
DiscML allows binary (1s/0s) analysis on data with more than two character states by converting all non-zero characters to 1s with simplify=TRUE. Bp Figure 1 Phylogenetic relationship of three Bacillaceae (B1, B2, B3) clades, on which the evolutionary rates of gene families are estimated using DiscML. A, a constant rate is estimated on each phylogeny; B, separate rates are estimated for external branches (μ 1 ) versus internal branches (μ 2 ) on each phylogeny. These three clades were studied in our previous study on gene presence, absence, and fragments [20]. Gene families are recategorized, with gene absence and fragments as character state 0, single-copy genes as 1, and gene families with two or more members as 2.

Results and discussion
DiscML was first tested using the gene family data on three Bacillaceae clades ( Figure 1A, Additional file 1 and [20]). In the previous study [20], we distinguished gene fragments from gene absence and gene presence. In this study, we eliminated the character state specific for gene fragments and re-categorized gene fragments as gene absence or character state 0, single-copy genes as character state 1, and gene families with two or more members 2 LnL 732 * * * 14 * * * 20 * * * μ 1 is for external branches, while μ 2 is for internal branches on each tree as illustrated in Figure 1B. * * * P < 0.001 (df=1), as 2 LnL approximately follows a χ 2 -distribution. http://www.  as 2 (Additional file 1), so that the application of BD models on these data is meaningful. It is worth to note that, though the number of character states is restricted to three here, DiscML is flexible and capable of analyzing a large number of character states. The performance of DiscML is found to be reliable. For instance, the ER+0 model with the option of simplify=TRUE in Table 1 is mathematically identical to the M 00 model in [20]. The optimization in [20] was achieved using the Nelder-Mead simplex method [25], while the optimization in Table 1 was achieved using the PORT routines [17]. Importantly, the DiscML estimates are identical to the previous estimates for all three clades. As expected, the parameter-rich models consistently outperformed the nested simplistic models (e.g., LnL BDARD > LnL BDISYM > LnL BDER; LnL ARD > LnL SYM > LnL ER). Consistent with previous studies [3,20,26], rate estimates in closely related clades tend to be higher than those in distantly related clades due to the transient nature of many acquired genes (Table 1). Tested on an Intel Core i7 (3.4 Ghz) 16 GB RAM Dell desktop, the computation using DiscML is fast (Table 2) (Table 2).
DiscML was developed to allow separate rates among branches since evolutionary rates can vary among lineages [27][28][29]. In the three Bacillaceae clades, we assigned separate rates between external branches (μ 1 ) and internal branches (μ 2 ) as illustrated in Figure 1B. Our results in Table 3 support the previous findings of higher gene turnover rates on external branches than those on internal branches [26,30].
It is often of interest for users to know the individual rate of each character site. Previously, we have shown that the mitochondrial intron in the 21S rRNA gene undergoes very rapid turnover in yeast [31]. In this study, we estimated the individual rates of all 17 mitochondrial introns on the yeast phylogeny ( Figure 2 and Additional file 1) based on the intron distribution pattern (Additional file 1). On the plot generated by DiscML using ind=TRUE (Figure 3), users can visually compare the individual rates of different introns. For instance, the introns at sites 7 and 8 have faster turnover rates than the 21S rRNA intron at site 17 ( Figure 3). The R commands used in the study are provided in Additional file 1.

Conclusion
We illustrated the versatility of DiscML on different types of data and analyses. With a great flexibility and  fast computational speed, we are confident that DiscML can be used in a variety of studies on different discrete characters.