Generation of the Exact Distribution and Simulation of Matched Nucleotide Sequences on a Phylogenetic Tree

Abstract

Nucleotide sequences are often generated by Monte Carlo simulations to address complex evolutionary or analytic questions but the simulations are rarely described in sufficient detail to allow the research to be replicated. Here we briefly review the Markov processes of substitution in a pair of matching (homologous) nucleotide sequences and then extend it to k matching nucleotide sequences. We describe calculation of the joint distribution of nucleotides of two matching sequences. Based on this distribution, we give a method for simulation of the divergence matrix for n sites using the multinomial distribution. This is then extended to the joint distribution for k nucleotide sequences and the corresponding 4 k divergence array, generalizing Felsenstein (Journal of Molecular Evolution, 17, 368–376, 1981), who considered stationary, homogeneous and reversible processes on trees. We give a second method to generate matched sequences that begins with a random ancestral sequence and applies a continuous Markov process to each nucleotide site as in Rambaut and Grassly (Computer Applications in the Biosciences, 13, 235–238, 1997); further, we relate this to an equivalent approach based on an embedded Markov chain. Finally, we describe an approximate method that was recently implemented in a program developed by Jermiin et al. (Applied Bioinformatics, 2, 159–163, 2003). The three methods presented here cater for different computational and mathematical limitations and are shown in an example to produce results close to those expected on theoretical grounds. All methods are implemented using functions in the S-plus or R languages.