A matrix-analytical sampling formula for time-homogeneous coalescent processes under the infinite sites mutation model

In this paper we develop a general framework for calculating the probability of a genetic sample under a time-homogeneous coalescent process and the infinite sites mutation model. The evolutionary model that we consider can be characterized as a two-step procedure: A coalescent process that describes the ancestral relatedness of the samples and a sprinkling of mutations in separate sites on the ancestral tree according to a Poisson process. The coalescent process is defined using multivariate phase-type theory. The requirements are a rate matrix that determines the transition rates between the ancestral states, an initial state probability vector, and a reward matrix that informs about the characteristics of the ancestral states. For example, the reward matrix could contain information about the number of singleton, doubleton or higher-order lineages in the ancestral states. We analyze the probability generating function for the evolutionary model as a function of the initial state probability vector, the transition rate matrix, the reward matrix, and the mutation rate. The matrix-analytical expression of the probability generating function allows us to develop a general method for calculating the probability of a population genetic data set. We demonstrate that the method is computationally attractive for a small number of mutations and provide a simple and easy-to-implement algorithm for determining the probability of a sample from the evolutionary model. The method is computationally stable and only involves a single inverse matrix operation, matrix multiplications and matrix additions. We provide comprehensive understanding of the procedure by detailed calculations and discussions of several elementary examples. These examples include different sample representations (labeled samples and the site frequency spectrum) and different demographic and genetic models (the structured coalescent and the Beta-coalescent). We apply the sampling formula to calculate probabilities of spectra for the Kingman coalescent and the Beta-coalescent. Even for a small number of samples and mutations we find that the probabilities for spectra vary in huge orders of magnitudes. We compare the probabilities of the spectra to the values of Tajima's D-statistics, and find that the D-statistic is a poor predictor for the probability of a spectrum. Finally, we investigate how the probabilities of the spectra vary with the parametrization of the Beta-coalescent.