2.1.1 Syntax.

In LPhy, each line’s syntax consists of a variable declaration on the left-hand side, a specification operator, and a generator on the right-hand side, with lines ending with a semicolon character. For instance:

In this example, variable b is specified by the normal generative distribution Normal() with two arguments: the mean mean (0.0) and the standard deviation sd (1.0).

The left-hand side declares the name of a variable or an array of variables (case sensitive). The right-hand side specifies the values of the variable or array. This can be a constant value, array of constant values, deterministic function or stochastic generative distribution. Deterministic functions and generative distributions are matched by method signatures from constructors of their corresponding Java class in the reference implementation. See the LPhy reference implementation manual available from the homepage https://linguaphylo.github.io/ for a complete list of functions and generative distributions. Arguments inside functions or generative distributions follow the convention: (argument name) = (value).

2.1.2 Specification operators.

An equal sign = is used to specify deterministic or constant values for variables, such as:

A tilde sign ∼ is a specification operator which denotes the relationship between a stochastic random variable and its generator. In a Bayesian context, this is the prior distribution of the random variable. Or when given observed data, this specifies the likelihood function. For example, this specifies a prior for the variable b:

The type of a variable is inferred from the return type of its generator and does not need to be declared. For arguments of functions or generative distributions, the types are defined in the LPhy reference implementation manual available from the homepage https://linguaphylo.github.io/.

2.1.3 Arrays.

Arrays can be defined using square brackets with elements delimited by comma separators. For sequences of consecutive integers, we allow a more compact notation using the colon ‘:’ to define a range. For example:

The variable c has an array with values (1, 2, 3, 4, 5). The variable d has an array with values 2 to 10.

2.1.4 Code blocks.

The data and model keywords are reserved to specify code blocks inside curly brackets. The data block is used to read in and store input data, which are used by the model. Within the data block, we can read in alignment data via the NEXUS or FASTA parsers, specify constant values, and store metadata about the dataset. The model block is used to define the models and parameters in a Bayesian phylogenetic analysis.

Example 1. An LPhy script defining a constant-size coalescent tree prior with log-normally distributed population sizes, a strict clock model, and a Jukes-Cantor model on 10 nucleotide sequences with 200 sites (base pairs).

Example 1 specifies a complete phylogenetic model using only five lines of code inside two blocks. The first block specifies ‘200’ nucleotide sites in ‘L’ and ten taxa named from 1 to 10 in ‘taxa’. Taxa can be declared as strings or as numbers. In this example, the taxa names are numbered. The second block declares stochastic nodes or random variables highlighted in bluishgreen, and their generative distributions highlighted in blue. Constant nodes with fixed values are shown in vermillion.

Fig 1 shows this model specification represented as a probabilistic graphical model. Stochastic nodes are shown as circles, deterministic functions are shown as diamonds, and constants are shown as squares.

2.1.5 Variable vectorization.

Named variables can be scalars or vectors. Any generator can be vectorized to produce a vector of IID random variables using the replicates keyword:

In the example above, κ is a random vector of three log-normally distributed IID values.

Vectorization can also be applied to a generative distribution that produces vectors. In this case, the output will be a matrix, as seen in the following example.

Here, π contains 3 vectors, where each vector represents nucleotide base frequencies. So the resultant matrix will be 3 x 4 (major dimension of 3 and a minor dimension of 4).

A second mechanism for vectorization can be used by passing a vector of elements instead of a single element as an input argument of a generator:

Here, κ and π are vectors, both with a major dimension of 3. These are passed as arguments into the hky deterministic function, which returns a vector containing three instantaneous rate matrices stored in Q.

2.1.6 Extension mechanism.

LPhy is designed to be modular and extensible to encourage integrative software development and adoption by method developers in the field. LPhy core and its extensions are built using Java 17 with long-term support (LTS) and Gradle 4. The LPhy extension mechanism is implemented using the Java Platform Module System (JPMS) and the Java Service Provider Interface (SPI). This modular extension framework allows developers greater flexibility in code development and software releases. New functionality such as generative distributions, data types, and deterministic functions can be developed within LPhy extension modules. Additionally, software releases can be done independently to core LPhy releases. An example of the LPhy extension mechanism is the Phylonco package [27] which is described in the results section.

2.1.7 Parametric distributions.

The LPhy reference implementation comes with a series of parametric distributions commonly used in evolutionary models, including uniform, normal, log-normal, gamma, exponential, and dirichlet distributions. Parametric distributions can be specified as generative distributions for model parameters by:

Each parametric distribution is characterized by its own parameters. In the example above, the extinction rate parameter μ is drawn from a log-normal distribution with mean -5 and standard deviation 1.25 in log space.

2.1.8 Tree models.

Tree models are used to generate phylogenetic trees, and are central components in phylogenetic simulation and analysis. We briefly describe some of the main tree models implemented in LPhy below. This includes coalescent models and birth-death models.

Serially sampled coalescent model

The simplest coalescent model LPhy implements is the constant-population size coalescent, which can be extended to generate serially sampled (heterochronous) data [28]:

The script above specifies a serially sampled constant-population size coalescent (a generative distribution) for tree ψ with four taxa, ‘a’, ‘b’, ‘c’, ‘d’ sampled at 0.0, 1.0, 2.0 and 3.0 time points, respectively. Here, sample time is defined as the age of a sample, where 0.0 represents the present moment.

Structured coalescent model.

The structured coalescent [29, 30] generalizes the constant-population size coalescent [25] by allowing multiple demes, each of which are characterized by a distinct population size. In the simplest case this population size does not change through time. Demes exchange individuals according to migration rates m specified in the off-diagonal elements of a migration matrix M, where the diagonal elements store the population sizes, θ, of each deme. For K demes, the population size parameter ‘theta’ is a K-tuple, and m is a (K² − K)-tuple.

In the example above, migrationMatrix is a deterministic function and StructuredCoalescent is a generative distribution. A stochastic node ‘g’ stores a gene tree sampled from a two-deme structured coalescent process.

Skyline coalescent model.

The skyline coalescent model [31] is a coalescent process that models changes in population sizes. This model is characterized having a constant population size for each coalescent interval, with instantaneous changes in population size at some coalescent events.

The following script specifies a Skyline coalescent model with 10 coalescent intervals (11 taxa), with four distinct population sizes.

Here, ‘g’ is a stochastic node in the PGM, with its value sampled from the SkylineCoalescent generative distribution. Ten coalescent intervals are defined through the ‘groupSizes’ argument: the first four coalescent intervals will be drawn assuming a ‘theta’ of 0.1, the next three intervals with ‘theta’ equal to 0.2, and so on.

Birth-death models.

Birth-death models are commonly used in macroevolution as sampling distributions for species trees. Models that parameterize the fossilization process can be especially useful, as they allow users to leverage fossil ages as data. When fossil morphological characters have also been scored, total-evidence dating can be carried out [12]. One such tree model is the serially sampled birth-death process [32], whose parameters ‘psi’ and ‘rho’ below represent the rate of sampling extinct and extant lineages, respectively:

Other tree models include the birth-death [33] and fossilized birth-death processes [34], as well as the Yule process [35].

2.1.9 Substitution models.

Substitution models consist of continuous-time Markov chains (CTMC) used to model the evolution of discrete characters, such as nucleotides and amino acid residues. LPhy implements a general formulation of a phylogenetic CTMC, known as the GTR model [36], under which several nested models can be specified. The first line below constructs the instantaneous rate matrix Q for an HKY model [37], which is then used to in PhyloCTMC, the generative distribution for the sequence alignment data D:

Other substitution models can be easily specified by passing different instantaneous transition rate matrices Q to PhyloCTMC, e.g., the matrix of the Jukes-Cantor model [24]:

For forward simulation PhyloCTMC is used as a generative distribution for a multiple sequence alignment, which is here represented by stochastic node ‘D’. When the model is employed for statistical inference, and data D is known, the PhyloCTMC represents the phylogenetic likelihood. More details are discussed in the Data clamping section.

2.1.10 Evolutionary clock models.

Molecular clocks are used to model the rate of evolutionary change and how it varies over time. The LPhy language supports strict clock [38, 39], local clock [40] and relaxed clock [41] models. Specifying a clock model is done by generating evolutionary rate values, one per phylogenetic tree branch, and then multiplying those rates by the length of the corresponding branch. The branches of the tree are measured in units of time, effectively scaling the tree to the number of expected substitutions per site.

The simplest clock model is the strict clock, where the evolutionary rate remains constant over the entire tree. Specifying a strict molecular clock can be done by specifying the ‘mu’ parameter in the PhyloCTMC distribution. The default value for the clock rate ‘mu’ is 1.0.

More realistic clock models like the uncorrelated relaxed clock model [41] assume the rate for each branch is drawn according to a parametric distribution. For example, a relaxed clock with rates drawn from a log-normal distribution can be constructed as follows:

Here, 30 rates are drawn independently from a log-normal distribution, and then each is assigned to one of the 30 branches of tree ψ.

2.1.11 Inference and data clamping.

In addition to simulation, LPhy allows users to use a specified model for inference. We have developed LPhyBEAST which provides support for inference with the BEAST2 engine.

To set up an inferential analysis with LPhy ‘data clamping’ is performed similar to the Rev language [10]. Data clamping involves associating an observed value with a random variable in the model, which is represented as a stochastic node in the probabilistic graphical model (PGM). By clamping data to a node, the user is informing the inference engine that the value of that particular variable is known and will be conditioned on for the purpose of inference.

In LPhy, data clamping can be accomplished using the ‘data block’. This block allows the user to specify the observed values of certain variables in the model, effectively clamping these variables to their observed values during inference. This is useful when working with real data, as it allows the user to incorporate the observed data into the analysis and improve the accuracy of the results.

In LPhy, data clamping can be achieved using the ‘data block’, for example:

Example 2. An LPhy script for phylodynamic analysis of a virus dataset containing Respiratory syncytial virus subgroup A (RSVA) genomic samples [42, 43].

In Example 2, we used a Respiratory syncytial virus subgroup A (RSVA) dataset [42, 43] containing 129 molecular sequences coding for the G protein collected between years 1956 and 2002. We use three partitions corresponding to the codon position, an HKY substitution model [37], coalescent tree prior [25] and a strict molecular clock with a log-normal prior on the mean clock rate. Within the data block we clamp the value of ‘codon’, a stochastic node that appears below inside the model block. This is achieved by specifying a data node of the same name (codon) in the data block. In this example the data is vectorized into three codon positions to allow different site models for the different codon positions.