BACKGROUND
monte carlo simulation of sequence evolution is routinely used in assessing the performance of phylogenetic inference methods , multiple sequence alignment algorithms  and ancestral reconstruction .
monte carlo simulation of sequence evolution is also crucially important in the testing of competing evolutionary hypotheses  <cit> , yet the effect of insertions and deletions  is often ignored since the necessary tools were not available.
several software packages for simulating basic sequence evolution under popular substitution models have been published in the last decade, for example sdse <cit> , seq-gen <cit>  and the evolver program from the paml package  <cit> .
more recently published software goes beyond the limitations of earlier simulation tools, allowing, for example, the simulation of indel events, sequence regions evolving under different models/parameters, the use of non-homogeneous models allowing for different parameters on different evolutionary lineages  and the flexible simulation of genomic features  <cit> .
the r language  <cit>  is the leading open-source environment for statistical computing and graphics, extensively used in bioinformatics data analysis.
its use for the analysis of phylogenetic and evolutionary data is aided by the "analysis of phylogenetics and evolution"  package  <cit>  and a small ecosystem of packages extending its capabilities  <cit> .
the simulation of the evolution of continuous characters is possible using ape and discrete characters can be evolved along a tree according to an arbitrary rate matrix using the phangorn <cit>  and geiger <cit>  packages.
however, there is no r package currently supporting the simulation of indel events and sequence evolution with site-specific rates, nonsynonymous/synonymous rate ratios or other advanced features available in other phylogenetics software.
allowing for heterogeneous evolution is a fundamental part of virtually all modern phylogenetic analyses  <cit>  and realistic simulation of indel events is indispensable when benchmarking the performance of multiple alignment methods.
previous software does not handle indels realistically, posing potential problems for the downstream analyses.
most programs assume a uniform distribution of indel events across the simulated biological sequences, despite the fact that those are likely to have regions evolving under different selective constraints  <cit> .
some tools try to address this problem by allowing for partitions evolving under different models/parameters.
however, the deletions are often not allowed to cross partition boundaries, which creates an unrealistic "edge effect".
the correlation between the selective constraints on indels and substitution  <cit>  is another aspect of sequence evolution which cannot be handled properly just by defining partitions.
here we present phylosim, an object-oriented framework enabling the realistic monte carlo simulation of sequence evolution.
phylosim significantly extends the range of realistic evolutionary patterns that can be simulated, and is freely extensible within the r environment.
implementation
the phylosim framework - written in pure r - builds on the ape package and aims to complement it.
it also uses the r.oo package  <cit> , which provides class-object-oriented facilities with references on top of the default function-object-oriented framework, and depends on the compoisson and ggplot <dig> packages.
the released packages are freely available under the gnu general public licence version  <dig> from cran  <cit>  and the package download page  <cit> .
the package sources are also available from the phylosim github repository  <cit> .
RESULTS
phylosim uses the gillespie algorithm  <cit>  as a unified framework for simulating substitutions and other events such as insertions and deletions .
sequence evolution along a branch is simulated in two steps, iterated repeatedly: sampling the time of occurrence of the next event and then modifying the sequence object according to a randomly selected event.
the rate of occurrence of the next event is equal to the sum of all possible event rates, while the event to be performed is selected with a probability proportional to its rate.
after performing the event, the set of possible events is updated.
these steps are repeated until the available time  is exhausted.
as in the case of previous software  <cit> , time is defined in terms of expected substitutions per site and the neutral rates of all other processes are specified relative to that.
selective constraints on different types of events  can be incorporated in a natural way in the framework described above by accepting/rejecting the selected event with a probability determined by some of its characteristics .
the key features offered by phylosim are the following:
i. simulation of the evolution of a set of discrete characters with arbitrary states evolving by a continuous-time markov process with an arbitrary rate matrix.
ii.
explicit implementations of the most popular nucleotide, amino acid and codon substitution models.
iii.
the possibility to simulate evolution by a combination of substitution processes with arbitrary rate matrices acting on the same site.
iv.
simulation under the popular models of among-sites rate variation, such as the gamma  and invariant sites plus gamma  models.
v. the possibility to simulate with arbitrarily complex patterns of among-sites rate variation by setting the site-specific rates according to any r expression.
vi.
simulation with one or more separate insertion and deletion processes acting on the sequences, each sampling indel lengths from an arbitrary discrete distribution or an r expression .
vii.
all the rate variation features listed above  can be readily applied to modify the rates whereby insertion and deletion processes initiate events at given sites.
viii.
simulation of the effects of spatially variable functional constraints by site- and process-specific insertion and deletion tolerance parameters, which determine the rejection probability of a proposed insertion or deletion ; rescaled deletion processes speed up simulation when deletions are strongly selected against .
field indel models allow for the fine-grained control of selective constraints on indels and, unlike the partition approach, do not suffer from "edge effect" artifacts.
ix.
the possibility of having different processes and site- and process-specific parameters for every site, which allows for an arbitrary number of partitions in the simulated data.
x. simulation of heterotachy and other cases of time-non-homogeneous evolution by allowing the user to set "node hook" functions altering sites' properties at internal nodes of the phylogeny.
xi.
full control over the properties of the inserted sequences, which makes it possible to easily extend phylosim with new kinds of insertion processes, .
the validity of the framework has been tested by simulating the evolution of nucleotide, amino acid and codon sequences of increasing length and estimating the value of model parameters and branch lengths from the resulting alignments using the paml package  <cit> .
the results are summarized in appendix a  along with the computing time needed for simulation and estimation.
implementation using r naturally affects the amount of computing time and memory needed for the simulations, but we believe that this is balanced out by the unparalleled versatility offered by the framework.
phylosim is provided with extensive documentation.
in addition, a package 'vignette'  gives a series of examples illustrating the simulation of successively more complex evolutionary scenarios, from very simple and familiar models through to complicated heterogeneous evolutionary dynamics not available with other software.
further details of the field deletion models
a natural way to incorporate deletions into the gillespie framework is to assign an individual rate to every possible deletion event.
modelling in this manner is extremely general but requires a lot of specification: not only individual sites' tolerance to deletion but also of how they interact with neighbouring sites.
instead we propose a more restricted "field model" of deletion that generalises previous simple approaches to allow the rate at which deletions occur to vary across the sequence but only requires one parameter per site - its deletion tolerance - to be specified.
under this model, deletions are proposed in same manner as other events, specifying a rate of occurrence and a distribution of lengths, and then accepted or rejected based on sites they propose to remove.
firstly consider only single-site deletions and let each site, i, in the sequence have an associated deletion tolerance parameter, di ∈  <cit> , representing the probability that it is actually deleted given that a deletion is proposed.
sites where di =  <dig> are deleted at the background rate, sites with di <  <dig> are deleted more slowly, and sites with di =  <dig> are never deleted.
for proposed deletions that span multiple sites, ℐ, each site is considered independently and the proposed deletion is accepted if and only if every site accepts it: the total probability of acceptance is therefore Πi∈ℐ di.
this scheme allows functionally important "undeletable" sites and regions to be modelled, as well as the phenomenon of deletion hotspots.
it is natural to think of the background rate of deletion as a neutral rate but this is not necessary and can lead to the gillespie algorithm becoming inefficient: for example, an extremely deletion intolerant sequence will reject almost all deletions proposed and so waste many steps.
instead we can rescale the process and the deletion tolerances  so that deletions are proposed at a rate equal to what would occur if the entire sequence had a deletion tolerance equal to its most tolerant site.
an example: annotating a simulated alignment using prank
simulating sequence evolution is crucial when benchmarking any method which relies on the heterogeneity of the evolutionary signal in multiple alignments .
as an example of a potential use of the phylosim package, we simulated the evolution of a genomic region containing a small gene with two exons , which could be a practical way to assess the sensitivity of the genomic structure model  <cit>  implemented in the prank phylogeny-aware multiple alignment tool  <cit> .
we simulated the evolution of the genomic region along a phylogenetic tree of nine mammal species .
for added realism, we included in the simulation features like fixed start codons and splice sites, and a substitution process acting on the three functionally equivalent stop codons .
the r script used for the simulation  can be found in the examples directory of the package source repository  <cit> .
we used the webprank server  <cit>  to align the simulated sequences, and regarded an alignment position to be annotated as coding if the reported posterior probability of any of the three coding states was greater than  <dig> .
we transferred back the annotation to the "true" simulated multiple alignment through the human sequence and compared it to the true structure of the simulated region .
we found that the exons inferred by prank show a good overlap with the true simulated exons.
CONCLUSIONS
with the features listed above, phylosim permits simulations encompassing a wide range of complexity , from those involving simple indel models similar to tkf <dig>  <cit>  to realistic simulations of protein sequences containing domains with distinct characteristics as well as of whole genomic regions harbouring coding sequences with intron-exon structures .
extensibility is the most prominent feature of the framework, its design making very simple the implementation of new processes embodying novel events  and the adaptation of the simulator to whatever is required.
availability of complex evolutionary processes in different simulation software.
additional details for less advanced software and simpler models are given by .
*see text for details.
availability and requirements
• project name: phylosim
• project home page: http://www.ebi.ac.uk/goldman-srv/phylosim
• project source repository: http://github.com/sbotond/phylosim
• operating system: os independent
• programming language: r
• required r packages: r.oo , ape , compoisson , ggplot <dig>
• license: gnu general public license version 3
• any restrictions to use by non-academics: none
authors' contributions
bs, tm and ng designed the framework.
bs implemented the framework and obtained the test results.
gj contributed the alignment and tree plotting methods.
bs drafted the manuscript, which was reviewed and approved by all authors.
supplementary material
additional file 1
the phylosim package vignette.
click here for file
additional file 2
appendix a.
click here for file
acknowledgements
bs was funded by an embo short-term fellowship and an embl interdisciplinary postdoc  under marie curie actions .
gj was funded by a gates cambridge scholarship and is a member of darwin college, university of cambridge.

 monte carlo simulation of sequence evolution is routinely used in assessing the performance of phylogenetic inference methods, multiple sequence alignment algorithms and ancestral reconstruction. monte carlo simulation of sequence evolution is also crucially important in the testing of competing evolutionary hypotheses, yet the effect of insertions and deletions is often ignored since the necessary tools were not available. several software packages for simulating basic sequence evolution under popular substitution models have been published in the last decade, for example sdse, seq-gen and the evolver program from the paml package. more recently published software goes beyond the limitations of earlier simulation tools, allowing, for example, the simulation of indel events, sequence regions evolving under different models/parameters, the use of non-homogeneous models allowing for different parameters on different evolutionary lineages and the flexible simulation of genomic features. the r language is the leading open-source environment for statistical computing and graphics, extensively used in bioinformatics data analysis. its use for the analysis of phylogenetic and evolutionary data is aided by the "analysis of phylogenetics and evolution" package and a small ecosystem of packages extending its capabilities. the simulation of the evolution of continuous characters is possible using ape and discrete characters can be evolved along a tree according to an arbitrary rate matrix using the phangorn and geiger packages. however, there is no r package currently supporting the simulation of indel events and sequence evolution with site-specific rates, nonsynonymous/synonymous rate ratios or other advanced features available in other phylogenetics software. allowing for heterogeneous evolution is a fundamental part of virtually all modern phylogenetic analyses and realistic simulation of indel events is indispensable when benchmarking the performance of multiple alignment methods. previous software does not handle indels realistically, posing potential problems for the downstream analyses. most programs assume a uniform distribution of indel events across the simulated biological sequences, despite the fact that those are likely to have regions evolving under different selective constraints. some tools try to address this problem by allowing for partitions evolving under different models/parameters. however, the deletions are often not allowed to cross partition boundaries, which creates an unrealistic "edge effect". the correlation between the selective constraints on indels and substitution is another aspect of sequence evolution which cannot be handled properly just by defining partitions. here we present phylosim, an object-oriented framework enabling the realistic monte carlo simulation of sequence evolution. phylosim significantly extends the range of realistic evolutionary patterns that can be simulated, and is freely extensible within the r environment. implementation the phylosim framework - written in pure r - builds on the ape package and aims to complement it. it also uses the r.oo package, which provides class-object-oriented facilities with references on top of the default function-object-oriented framework, and depends on the compoisson and ggplot packages. the released packages are freely available under the gnu general public licence version from cran and the package download page. the package sources are also available from the phylosim github repository.  phylosim uses the gillespie algorithm as a unified framework for simulating substitutions and other events such as insertions and deletions. sequence evolution along a branch is simulated in two steps, iterated repeatedly: sampling the time of occurrence of the next event and then modifying the sequence object according to a randomly selected event. the rate of occurrence of the next event is equal to the sum of all possible event rates, while the event to be performed is selected with a probability proportional to its rate. after performing the event, the set of possible events is updated. these steps are repeated until the available time is exhausted. as in the case of previous software, time is defined in terms of expected substitutions per site and the neutral rates of all other processes are specified relative to that. selective constraints on different types of events can be incorporated in a natural way in the framework described above by accepting/rejecting the selected event with a probability determined by some of its characteristics. the key features offered by phylosim are the following: i. simulation of the evolution of a set of discrete characters with arbitrary states evolving by a continuous-time markov process with an arbitrary rate matrix. ii. explicit implementations of the most popular nucleotide, amino acid and codon substitution models. iii. the possibility to simulate evolution by a combination of substitution processes with arbitrary rate matrices acting on the same site. iv. simulation under the popular models of among-sites rate variation, such as the gamma and invariant sites plus gamma models. v. the possibility to simulate with arbitrarily complex patterns of among-sites rate variation by setting the site-specific rates according to any r expression. vi. simulation with one or more separate insertion and deletion processes acting on the sequences, each sampling indel lengths from an arbitrary discrete distribution or an r expression. vii. all the rate variation features listed above can be readily applied to modify the rates whereby insertion and deletion processes initiate events at given sites. viii. simulation of the effects of spatially variable functional constraints by site- and process-specific insertion and deletion tolerance parameters, which determine the rejection probability of a proposed insertion or deletion; rescaled deletion processes speed up simulation when deletions are strongly selected against. field indel models allow for the fine-grained control of selective constraints on indels and, unlike the partition approach, do not suffer from "edge effect" artifacts. ix. the possibility of having different processes and site- and process-specific parameters for every site, which allows for an arbitrary number of partitions in the simulated data. x. simulation of heterotachy and other cases of time-non-homogeneous evolution by allowing the user to set "node hook" functions altering sites' properties at internal nodes of the phylogeny. xi. full control over the properties of the inserted sequences, which makes it possible to easily extend phylosim with new kinds of insertion processes,. the validity of the framework has been tested by simulating the evolution of nucleotide, amino acid and codon sequences of increasing length and estimating the value of model parameters and branch lengths from the resulting alignments using the paml package. the results are summarized in appendix a along with the computing time needed for simulation and estimation. implementation using r naturally affects the amount of computing time and memory needed for the simulations, but we believe that this is balanced out by the unparalleled versatility offered by the framework. phylosim is provided with extensive documentation. in addition, a package'vignette' gives a series of examples illustrating the simulation of successively more complex evolutionary scenarios, from very simple and familiar models through to complicated heterogeneous evolutionary dynamics not available with other software. further details of the field deletion models a natural way to incorporate deletions into the gillespie framework is to assign an individual rate to every possible deletion event. modelling in this manner is extremely general but requires a lot of specification: not only individual sites' tolerance to deletion but also of how they interact with neighbouring sites. instead we propose a more restricted "field model" of deletion that generalises previous simple approaches to allow the rate at which deletions occur to vary across the sequence but only requires one parameter per site - its deletion tolerance - to be specified. under this model, deletions are proposed in same manner as other events, specifying a rate of occurrence and a distribution of lengths, and then accepted or rejected based on sites they propose to remove. firstly consider only single-site deletions and let each site, i, in the sequence have an associated deletion tolerance parameter, di ∈, representing the probability that it is actually deleted given that a deletion is proposed. sites where di = are deleted at the background rate, sites with di < are deleted more slowly, and sites with di = are never deleted. for proposed deletions that span multiple sites, ℐ, each site is considered independently and the proposed deletion is accepted if and only if every site accepts it: the total probability of acceptance is therefore Πi∈ℐ di. this scheme allows functionally important "undeletable" sites and regions to be modelled, as well as the phenomenon of deletion hotspots. it is natural to think of the background rate of deletion as a neutral rate but this is not necessary and can lead to the gillespie algorithm becoming inefficient: for example, an extremely deletion intolerant sequence will reject almost all deletions proposed and so waste many steps. instead we can rescale the process and the deletion tolerances so that deletions are proposed at a rate equal to what would occur if the entire sequence had a deletion tolerance equal to its most tolerant site. an example: annotating a simulated alignment using prank simulating sequence evolution is crucial when benchmarking any method which relies on the heterogeneity of the evolutionary signal in multiple alignments. as an example of a potential use of the phylosim package, we simulated the evolution of a genomic region containing a small gene with two exons, which could be a practical way to assess the sensitivity of the genomic structure model implemented in the prank phylogeny-aware multiple alignment tool. we simulated the evolution of the genomic region along a phylogenetic tree of nine mammal species. for added realism, we included in the simulation features like fixed start codons and splice sites, and a substitution process acting on the three functionally equivalent stop codons. the r script used for the simulation can be found in the examples directory of the package source repository. we used the webprank server to align the simulated sequences, and regarded an alignment position to be annotated as coding if the reported posterior probability of any of the three coding states was greater than. we transferred back the annotation to the "true" simulated multiple alignment through the human sequence and compared it to the true structure of the simulated region. we found that the exons inferred by prank show a good overlap with the true simulated exons.  with the features listed above, phylosim permits simulations encompassing a wide range of complexity, from those involving simple indel models similar to tkf to realistic simulations of protein sequences containing domains with distinct characteristics as well as of whole genomic regions harbouring coding sequences with intron-exon structures. extensibility is the most prominent feature of the framework, its design making very simple the implementation of new processes embodying novel events and the adaptation of the simulator to whatever is required. availability of complex evolutionary processes in different simulation software. additional details for less advanced software and simpler models are given by. *see text for details. availability and requirements • project name: phylosim • project home page: http://www.ebi.ac.uk/goldman-srv/phylosim • project source repository: http://github.com/sbotond/phylosim • operating system: os independent • programming language: r • required r packages: r.oo, ape, compoisson, ggplot • license: gnu general public license version 3 • any restrictions to use by non-academics: none authors' contributions bs, tm and ng designed the framework. bs implemented the framework and obtained the test results. gj contributed the alignment and tree plotting methods. bs drafted the manuscript, which was reviewed and approved by all authors. supplementary material additional file 1 the phylosim package vignette. click here for file additional file 2 appendix a. click here for file acknowledgements bs was funded by an embo short-term fellowship and an embl interdisciplinary postdoc under marie curie actions. gj was funded by a gates cambridge scholarship and is a member of darwin college, university of cambridge.
the monte carlo simulation of sequence evolution is routinely used to assess the performance of phylogenetic inference methods and sequence alignment algorithms. phylosim is an extensible framework for the monte carlo simulation of sequence evolution, written in r, using the gillespie algorithm to integrate the actions of many concurrent processes such as substitutions, insertions and deletions. close integration with r and the wide range of features implemented offer unmatched flexibility, making it possible to simulate sequence evolution under a wide range of realistic settings. uniquely among sequence simulation tools, phylosim can simulate arbitrarily complex patterns of rate variation and multiple indel processes, and allows for the incorporation of selective constraints on indel events. progress in the field of molecular evolution fuels the need for more realistic and hence more complex simulations, adapted to particular situations, yet current software makes unreasonable assumptions such as homogeneous substitution dynamics or a uniform distribution of indels across the simulated sequences. user-defined complex patterns of mutation and selection can be easily integrated into simulations, allowing phylosim to be adapted to specific needs.
 the monte carlo simulation of sequence evolution is routinely used to assess the performance of phylogenetic inference methods and sequence alignment algorithms. progress in the field of molecular evolution fuels the need for more realistic and hence more complex simulations, adapted to particular situations, yet current software makes unreasonable assumptions such as homogeneous substitution dynamics or a uniform distribution of indels across the simulated sequences. this calls for an extensible simulation framework written in a high-level functional language, offering new functionality and making it easy to incorporate further complexity.  phylosim is an extensible framework for the monte carlo simulation of sequence evolution, written in r, using the gillespie algorithm to integrate the actions of many concurrent processes such as substitutions, insertions and deletions. uniquely among sequence simulation tools, phylosim can simulate arbitrarily complex patterns of rate variation and multiple indel processes, and allows for the incorporation of selective constraints on indel events. user-defined complex patterns of mutation and selection can be easily integrated into simulations, allowing phylosim to be adapted to specific needs.  close integration with r and the wide range of features implemented offer unmatched flexibility, making it possible to simulate sequence evolution under a wide range of realistic settings. we believe that phylosim will be useful to future studies involving simulated alignments. 