## Abstract

Most studies of mutation rates implicitly assume that they remain constant throughout development of the germline. However, researchers recently used a novel statistical framework to reveal that mutation rates differ dramatically during sperm development in *Drosophila melanogaster*. Here a general framework is described for the inference of germline mutation patterns, generated from either mutation screening experiments or DNA sequence polymorphism data, that enables analysis of more than two mutations per family. The inference is made more rigorous and flexible by providing a better approximation of the probabilities of patterns of mutations and an improved coalescent algorithm within a single host with realistic assumptions. The properties of the inference framework, both the estimation and the hypothesis testing, were investigated by simulation. The refined inference framework is shown to provide (1) nearly unbiased maximum-likelihood estimates of mutation rates and (2) robust hypothesis testing using the standard asymptotic distribution of the likelihood-ratio tests. It is readily applicable to data sets in which multiple mutations in the same family are common.

SPERM and eggs experience many divisions after the fertilized egg, and mutations may occur each time a cell divides. Little is known about the patterns of mutations during development of the germline cell lineage. This is partly due to the scarcity of appropriate experimental data and to lack of proper statistical methods for analyzing such data. Recently, Gao *et al.* (2011) reported that mutation rates differ dramatically during germline development in *Drosophila*, with the rate for the first cell division the highest. But the method developed by Gao *et al.* (2011) is limited to handling only families with at most two mutations each. Also, their conclusions relied heavily on hypothesis testing and the statistical properties of the likelihood-ratio test they used are not known under such circumstances. Furthermore, their coalescent algorithm is too simplistic, not taking into consideration the details of spermatogenesis. Since the ability to make inferences about mutation rates at the level of single-cell division would be a significant step forward, it is desirable to make the inference rigorous and applicable for analysis of data in which more than two mutations per family are common.

Knowledge of development of the germline lineage is essential for inferring mutation rates. For *Drosophila melanogaster* males, each sperm from a young adult has experienced ≥36 divisions. The first 14 divisions occur in the cleavage stage characterized by fast cell divisions; the last 5 occur during spermatogenesis; those in between occur during gastrulation and organogenesis when the germline stem cells (GSC) divide asymmetrically. For *Drosophila*, it is well known (Drost and Lee 1995, 1998; Gilbert 2003) that (1) after the 8th cell division, ∼4–6 cells become the primordial germ cells (PGC); (2) after the 12th division, the number of PGCs ranges from 23 to 52; (3) after the 14th division there are 5–6 PGCs in each gonad; and (4) from the 15th division to just before spermatogenesis the number of PGCs remains more or less constant. After the 31st division one of the two daughter cells of each stem cell remains as a stem cell; the other one differentiates into 64 sperm. Thus, if sperm are sampled after the 36th division, all have experienced exactly 36 divisions, but if they are sampled after, for example, the 38th division, some would have experienced 36 divisions, some 37, and some 38 divisions. The algorithm developed by Gao *et al.* (2011), using the principle of coalescence (Kingman 1982; Ewens 2004), does not account for these differences.

To develop a thorough understanding of mutational patterns during germline development requires obtaining an estimate of mutation rate for each cell division and testing various hypotheses about mutation rates. This article describes further development of the inference framework, which overcomes previous shortcomings, investigates its statistical properties, improves the coalescent algorithm, and reanalyzes the published data. The improved inference framework has the advantage of being adaptable to analyzing mutation patterns in nucleotide polymorphism data, such as generated by next-generation DNA sequencing.

## The Theory

### Definitions and notations

Consider a sample of families, each consisting of a number of offspring (or sperm) from the same father. For each family with *n* sampled offspring, a *mutation pattern* is observed and represented by 〈*i*_{1}, *i*_{2}, … , *i _{k}*〉 such that each element represents an identified mutation, its value equal to the number of mutants for the mutation, where

*k*is the number of mutations. For example, 〈1〉 represents a mutation in a family that leads to only one mutant; 〈3, 2〉 represents two mutations, where one leads to three mutants and another to two mutants. Also we use 〈〉 to represent the case in which no mutation is observed.

Sequencing the same region of the genome among all sampled offspring in the same family will yield a mutation pattern. Alternatively, such information can be obtained from traditional experiments, particularly for some model organisms. For *Drosophila*, multigeneration mutation screening has been well developed (Muller 1928; Woodruff *et al.* 1984, 1996; Ashburner 1989; Greenspan 1997) and one such system was used by Gao *et al.* (2011) for the purpose of identifying recessive lethal or nearly lethal mutations. More experimental detail can be found in Gao *et al.* (2011) and for the main purpose of this article, it should suffice to outline the structure of information and the type of mutation being examined. The *recessive lethality d* of a recessive mutation is defined as one minus the maximal percentage of the homozygote (among all survival offspring) for that mutation. The data reported by Gao *et al.* (2011) correspond to those mutations with recessive lethality equal to 99%, that is, no more than 1% of survival offspring are z/z homozygote. Once a cell acquires a mutation with recessive lethality *d*, further mutation(s) is much more likely to increase lethality than to reverse it. Consequently recessive lethal mutations have a masking effect such that only the earliest one is identifiable. Figure 1 shows examples of how mutations in a genealogy lead to different mutation patterns.

Suppose branches of a sample genealogy are labeled by integers. For branch *i*, define Ω(*i*) as the set of branches consisting of the branch *i* and all its descendant branches, which are referred to as the subtree of the branch *i*. For the genealogy shown in Figure 2, for example, Ω(2) = {2, 5, 6} and Ω(3) = {3, 4, 7, 9, 8, 9}. Therefore, two mutations, respectively on branches *i* and *j*, are both observable if and only if Ω(*i*) ∩ Ω(*i*) = Ø.

Suppose the germ cell divisions from a fertilized egg to sperm are divided into *I* intervals. Let [*i*, *j*] represent the interval from the *i*th to the *j*th cell divisions. Suppose the mutation rate per cell division for the *l*th interval is *u _{l}* and define

**u**= (

*u*

_{1}, … ,

*u*)

_{I}*. For a given sample of sperm, there is a genealogy connecting them to the fertilized egg. Suppose each branch in the genealogy is identified by an unique integer (how branches are numbered is immaterial). Define for the*

^{T}*i*th branch,

*b*as the number of divisions it contains from the

_{ij}*j*th interval and

**b**

*as a vector with elements*

_{i}*b*,

_{ij}*j*= 1, … ,

*I*. That is,

**b**

*= (*

_{i}*b*

_{i}_{1}, … ,

*b*)

_{iI}*. Define*

^{T}*φ*(

*i*) as the

*size*of the

*i*th branch,

*i.e.*, the number of descendants of the branch that are observed in the sample, and (1)Then

*a*is the total number of cell divisions from the

_{kj}*j*th interval that are of size

*k*and the

*k*th element,

*t*, of

_{k}**t**is the total number of cell divisions from the

*k*th interval. For branch

*i*, let

**w**

*be the sum of lengths of all the branches in Ω(*

_{i}*i*), excluding the branch

*i*itself. That is, . Figure 2 illustrates the aforementioned quantities in a genealogy of five alleles taken after the fifth division. It follows that ,

**a**

_{2}=

**b**

_{2}+

**b**

_{4}= (1, 3, 0)

*,*

^{T}**a**

_{3}=

**b**

_{3}= (1, 1, 0)

*,*

^{T}**a**

_{4}= 0,

**a**

_{5}=

**b**

_{1}= (1, 0, 0), and an example of

**w**is that

**w**

_{4}=

**b**

_{7}+

**b**

_{8}and

**w**

_{3}=

**b**

_{4}+

**b**

_{7}+

**b**

_{8}+

**b**

_{9}.

In addition to **a** and **t**, we will encounter other quantities that are functions of **b*** _{i}*,

*i*= 1, … , which will be defined as they are introduced. Each of these quantities has a value for a given genealogy, and often we need to evaluate its expectation (mean) over all genealogies. We use a bar over the variable to represent its expectation. For example,

### Probability of a mutation pattern

Assume that the number of mutations in a branch of the sample genealogy *g* is a Poisson variable. Then the probability of no mutation in a family is equal to . Since a single mutation leading to an observed pattern 〈*i*〉 must occur on a branch of size *i*, it follows that (3)where the summation is taken over all the branches of size *i*. In the summation, the first term is the probability that there is no mutation outside the subtree of branch *k* and the second term is the probability there is at least one mutation in branch *k*. This is because any mutation in the subtree will be masked by the mutation in branch *k* and thus not observable. In general, we have for a mutation pattern 〈*i*_{1}, … , *i _{l}*〉 that(4)where and

*(*

_{g}*i*

_{1}, … ,

*i*) is the collection of the branch sets of genealogy

_{l}*g*on which mutations can lead to the observed mutational pattern. That is,

Since sample genealogy is generally unobservable, one needs to consider all the possible sample genealogies from which the given mutational pattern can be generated, which leads to the general unconditional probability of the mutational pattern 〈*i*_{1}, … , *i _{l}*〉 as(6)This formula provides the basis for the proposed inferences and detailed analysis of

*Drosophila*data. For any given mutation pattern 〈

*i*

_{1}, … ,

*i*〉 and

_{l}**u**, the probability can be evaluated as the average of

*Pr*(〈

*i*

_{1}, … ,

*i*〉|

_{l}*g*) [which is given by (4)] over a reasonably large set of simulated sample genealogies. However, it is generally not efficient and often impractical to use the above formula directly if hundreds or even thousands of different

**u**need to be evaluated.

## Approximation to the Probability of a Mutation Pattern

Since Equation 6 is computationally expansive to use in general, accurate and yet-fast approximations to the probabilities of various mutation patterns are important and often necessary for large-scale data analysis. Gao *et al.* (2011) found approximations to the probabilities for up to two mutations in a family, using the Taylor expansion. For example, . Define . Then the probabilities *p*_{0} = *Pr*(〈〉), *p _{i}* =

*Pr*(〈

*i*〉), and

*p*=

_{ij}*Pr*(〈

*i*,

*j*〉) can be approximated (Gao

*et al.*2011) by (7) (8) (9)where

*δ*

_{i}_{−}

*= 1 if*

_{j}*i*=

*j*and 0 otherwise. This method of approximation is referred to as the approximation by Taylor expansion (ATE). Although these approximations cover up to two mutations per family, in principle a reasonably accurate approximation to the probability of any given mutation pattern can be obtained if a sufficient number of Taylor expansion terms are included. With increasing mutation rate, the number of required terms for each case will also increase, and due to the need to estimate a large number of coefficients in higher-order terms, their computations make the ATE inefficient.

Since typically in Equation 6, . Furthermore one can simplify the expression by replacing *w* for each combination of branches by its average value and arrive at (10) (11)where *w*(〈*i*_{1}, … , *i _{l}*〉) is defined as the average value of over the set

*(*

_{g}*i*

_{1}, … ,

*i*) and (12)where , which leads to an approximation of Equation 6 as(13)A further simplification and approximation can be obtained by moving the integration inward and replacing each quantity by its integral (that is, its expectation). This leads to the approximation (14)where is the mean value of

_{l}**w**(〈

*i*

_{1}, …,

*i*〉) over all genealogies, and is the mean of

_{l}*S*over all genealogies, which can be computed as (15)where is the mean of over all possible genealogies. This new approach is referred to as the approximation by inward integration (AII).

Let *S*(〈〉, *u*) = 1 and *w*({}) = 0. Then Equation 14 is applicable to any mutation pattern. The computation of Equation 14 is quite manageable now. In particular, for up to two mutations, we have (16) (17)where . Therefore, the corresponding new approximations up to two mutations are (18) (19) (20)where is the mean **B*** _{ij}* over all genealogies. Note that

**B**

*is not the same as*

_{ij}**A**

*due to constraints on the pair of branches that are compatible with the observed pattern. Gao*

_{ij}*et al.*(2011) recognized the masking effect of mutations and estimated

**A**

*by*

_{ij}**B**

*. We show in a later section that when three or more mutations in a family are rare, then both the ATE and the AII give excellent approximations to the true probabilities. With an increasing number of families with more than two mutations, we find that the new approach provides a more accurate approximation to Equation 6 than those by Gao*

_{ij}*et al.*(2011).

## The Likelihood Inference

Suppose there are in total *m* different mutation patterns in the data set, *c*_{1}, … , *c _{n}*, and

*n*is the occurrence of pattern

_{i}*c*. Then, the likelihood of the data is (21)where

_{i}*Pr*(

*c*) is the probability of pattern

_{i}*c*. Based on the new scheme for estimating the probabilities of each pattern, the maximum-likelihood estimates, , of

_{i}**u**can be derived from ln(

*L*), which is(22)The asymptotic covariance of the estimates can also be obtained as the inverse of matrix , whose computation is described in the

*Appendix*. Let

**r**

*= (*

^{T}*r*

_{1}, … ,

*r*), where

_{I}*r*is the number of cell divisions in the

_{k}*k*th interval. Then per generation mutation rate

*u*can be estimated as(23)The variance of this estimate is . Suppose the total number of mutant lines in the experiment is

*M*and the total number of lines screened is

*N*. Then an alternative estimate of

*u*is , which is unbiased regardless of whether mutation rates during development are the same (Fu and Huai 2003). A hypothesis can be tested through the likelihood-ratio test. For example, for testing the null hypothesis (H

_{0}) that mutation rates at different cell divisions are all equal, against the alternative hypothesis H

_{1}that rates have no constraint, the test statistic (24)follows asymptotically the

*χ*

^{2}-distribution with

*I*− 1 d.f.

## Cell Propagation and Simulation of Cell Genealogy

A discrete generation model is used for the propagation of cells in the germline lineage. We introduce two alternative modes of cell propagation.

Let *N*(*i*) be the size of the *i*th population, which can be divided into two groups, one [size *N*_{0}(*i*)] without a sister cell and one [size *N*_{1}(*i*)] with one sister cell [*N*(*i*) = *N*_{0}(*i*) + *N*_{1}(*i*)]. The first mode of propagation assumes that for each cell in the *i*th population, the probability of having *k*(*k* = 0, 1, 2) daughter cell(s) in the *i*th population is *p _{k}*.

This mode of cell propagation is fully determined when the values of *p _{i}* are specified.

*p*

_{2}= 1 corresponds to the case in which each cell yields two daughter cells, which is considered to be the default situation. Another special case is that every cell produces at least one daughter cell, which corresponds to

*p*

_{0}= 0. For two randomly selected cells from the

*i*th population, the probability that they will coalesce in the

*i*− 1th population is(25)That is, two cells will coalesce if and only if the first cell selected has a sister cell [with probability

*N*

_{1}(

*i*)/

*N*(

*i*)] and the second cell selected is its sister cell [with probability 1/(

*N*(

*i*− 1) − 1]. When there are multiple pairs of cells being considered, multiple coalescence can occur, which is usually not allowed in the conventional coalescent theory. The exact probabilities of any particular pattern of coalescence (for example, two pairs of coalescence, five pairs of coalescence, etc.) can be given analytically although they are not necessary for our purpose. What is critical is a proper algorithm to simulate this process as is discussed later in this section.

An alternative mode of cell propagation is as follows. Assume that each cell in the (*i* − 1)th population divides to yield two daughter cells and the cells in the *i*th population are a random sample (without replacement) from these 2*N*(*i* − 1) daughter cells. This mode of cell propagation is recognized when a range condition, such as *N*(*i*) ∈ [*a*, *b*], is specified. In such a case, *N*(*i*) is assumed to be a random integer in the given range [*a*, *b*]. Then the probability that two randomly selected cells from the *i*th population will coalesce in the (*i* − 1)th population is (26)which occurs only if the second cell selected is the sister cell of the first one. Again the probability of multiple coalescence can be derived. However, the sampling process will also yield *N*(*i*) and *N*_{1}(*i*); thus the coalescent probability is also given by Equation 25.

### Simulation algorithms

A forward–backward two-step algorithm was used in Gao *et al.* (2011) and will continue to be used here. The first (forward) step is to simulate a history of the population sizes and the second (backward) step is to simulate the genealogy given the history of the populations sizes as follows.

#### Forward algorithm: Simulation of cell population dynamics:

Given the value of *N*(*i* − 1), the value of *N*(*i*) is simulated according to the transition mode. Meanwhile, the values of *N _{k}*(

*i*), (

*k*= 1, 2) are recorded.

#### Backward algorithm: Simulation of cell genealogy:

Given a collection of *n* cells from the *i*th population:

Create an array of

*N*(*i*) integers as follows:*N*_{1}(*i*) integers from 1 to*N*_{1}(1) and two copies of each integer from*N*_{1}(*i*) + 1 to*N*_{1}(*i*) +*N*_{2}(*i*).Take a random sample of size

*n*from the above array. If two integers in the sample are the same, a coalescent event occurs.Update the collection of cells and repeat steps 1 and 2 until the 0th population (the zygote) is reached.

In the forward step, the cell lineage splits into two subpopulations after the 14th division and enters the stem cell lineage, which is specified by mode 1 with *p*_{1} = 1 − *e*, *p*_{2} = 2 for small values of *e*. After the 31st division, each of the differentiated cells from stem cells goes into spermatogenesis, resulting in 64 sperm. After the 37th cell division, the sperm population consists of sperm that are derived from differentiated cells that have experienced different numbers of divisions. Therefore, the number of cell divisions for each of the sperm in a random sample can be different.

## Numerical Results

To investigate statistical properties of the inference framework, the 36 cell divisions are divided into four intervals: [1, 3], [4, 14], [15, 31], and [32, 36], representing, respectively, the early cleavage, late cleavage, the stem-cell stage, and the spermatogenesis stage. Situations with maximal cell divisions >36 are also considered, so that different sperm in a sample might have experienced different numbers of divisions (see Figure 3). In such situations, the meaning of the last two intervals needs to be modified. For example, if a maximum of 38 divisions is allowed, then the last interval corresponds to the last 5 cell divisions, which for some lineages are from the 34th to the 38th division, while for some they are from the 32nd to the 36th division; and the second interval thus includes the 15th division in whatever is not included in the last interval.

### Simulation of sample genealogy

As pointed out earlier, the process of simulating a sample genealogy is similar to that in Gao *et al.* (2011), with the exception of the gametogenesis stage. The process consists of forward and backward steps. The former is guided by a number of constraints about the population sizes mentioned in the Introduction. For example, after the 8th division, there are 256 cells from which only 4–7 cells become PGCs. Table 1 lists all the used constraints for population sizes during germline development. After the 31st division in the forward process, each of the differentiated cells will go into gametogenesis, which progresses through 5 additional divisions to produce 64 sperm. This aspect of the development is now explicitly modeled.

The backward process of the simulation is the same as that in Gao *et al.* (2011) except that the efficiency of the program has been improved. The resulting population dynamics with relation to the sample genealogy are illustrated in Figure 3. It is important to simulate a large number of genealogies from which the values of various coefficients in the inference framework can be obtained. To deal with up to four mutations in a family, we found that in general 250,000 genealogies are sufficient.

### Accuracy of the approximations to the probabilities of mutation patterns

Since the expected numbers of occurrences and the difference in the expected numbers of occurrences are critical to the statistical inference, we use the following index to measure the accuracy of the approximations, (27)where *P _{i}* is the exact probability for a mutation pattern with

*i*mutations, is its approximation, and

*N*is the number of families, which is set to 8625. The probabilities were estimated with 2 million simulated genealogies. Table 2 gives the results for several mutation rates. In all cases the AII is better than the ATE; the AII performs well for a wide range of mutation rates, including a very large mutation rate. The ATE appears to be sufficient for a mutation rate up to ∼

*u*× 10

^{−4}, but with a rate closer to 10 × 10

^{−4}, its errors becomes too large, in addition to not being able to handle more than two mutations. We focus on further studying the statistical properties of the AII because of its obvious superiority.

### Maximum-likelihood estimate of *u*

One major outcome of the inference is the maximum-likelihood estimate of the mutation rate **u**; thus it is important to understand the properties of the estimates, which were carried out using simulations. Table 3 shows the means and standard deviations of the maximum-likelihood estimates of *u*_{1}, *u*_{2}, *u*_{3}, and *u*_{4} for several cases. The results show that the maximum-likelihood estimates are slightly biased but the bias decreases with increase in family number, which is expected from the well-known properties of the maximum-likelihood method. The standard errors of estimating *u*_{1}, … , *u*_{4} differ from each other, with those for *u*_{3} and *u*_{4} being the smallest and that for *u*_{1} the largest. This pattern agrees with the fact that there are many more mutations that result in a smaller mutant size, most of which likely occurred during the third and fourth time intervals. As a result, there are more observations from these two intervals that lead to more accurate estimates.

Figure 4 shows the distributions of estimates of mutation rates corresponding to the first row of Table 3. Two obvious features from these distributions are as follows. The first is that with increased family number, each distribution becomes more concentrated around the true mutation rate. The second is that judging from the spread of the distributions, the quality of estimations for *u*_{3} and *u*_{4} is better than that for *u*_{1} and *u*_{2}. Among the four, the quality of estimating *u*_{1} is the poorest. These features agree well with the patterns of standard deviations in Table 3.

### Likelihood-ratio test

Being able to obtain maximum-likelihood estimates under different assumptions also allows us to examine the distribution of the likelihood-ratio test. Various hypotheses about the pattern of mutation rates can be tested, as reported in Gao *et al.* (2011); however, the following four are representative, one for each value of the degrees of freedom: H_{0}, rates are constant; H_{1}, the last three are the same; H_{2}, the first two are the same; and H_{3}, no constraint.

Let *L _{ij}* be the log-likelihood ratio statistics between the

*i*th and

*j*th hypotheses. When H

_{0}is true, it is expected that

*L*

_{0,1},

*L*

_{0,2}, and

*L*

_{0,3}follow asymptotically

*χ*

^{2}-distributions with 1, 2, and 3 d.f., respectively. In the simulations, constant mutation rates are used and for each simulated sample, the maximum likelihood under each hypothesis is found, which leads to the likelihood-ratio statistics. Table 4 shows the upper-tail critical values for these three statistics. Comparing these critical values with the critical values of

*χ*with 1, 2, and 3 d.f., respectively, indicates that the these critical values agree reasonably well with the asymptotic values with sample sizes as small as 500. The distributions of these statistics are given in Figure 5 for two different sample sizes, which shows the overall excellent agreement of the empirical density with asymptotic ones.

### Reanalysis of the data

The data being reanalyzed here consist of those presented in Table 1 of Gao *et al.* (2011) and 7 additional families, 3 of which have three mutations and 1 of which has four mutations, giving thus a total of 8,625 families. For convenience of comparison, we used the same division of intervals: [1, 1], [2, 2], [3, 14], [15, 31], and [32, 36]. Table 5 shows the maximum-likelihood estimates using both the ATE and the AII (for the sake of space, only the results for four of the eight hypotheses considered in Gao *et al.* 2011 are given), while Table 6 gives the results of the likelihood-ratio tests. Comparing the entries of the ATE in these tables to those in Tables 4 and 5 of Gao *et al.* (2011), one can see that the differences are minimal. Furthermore, comparing the estimates by the ATE to those by the AII shows that the differences are also minor in almost all cases. Therefore, the improved method does not change the conclusions made previously. These analyses also included the case in which 38 cell divisions were assumed. In such situations, the patterns of the likelihood-ratio tests (Table 6) suggest that the mutation rates for the second, third, and fourth intervals may also be different, although the evidence is only marginal.

While it is comforting that the reanalysis reinforces the conclusions made earlier, this should not be regarded as the AII lacking importance. When the number of families with more than two mutations increases, one can expect to see increasing differences and a more rigorous new method than the ATE will become necessary. To illustrate, we simulated sets of 8625 families with four intervals of cells [1, 1], [2, 14], [15, 33], and [34, 38], using two sets of mutation rates, one being equal rates for all the intervals and the other being one that produces mutational patterns resembling those from the experiment, which will be subjected to detailed analysis elsewhere. Table 7 shows the comparison of the two methods from which it is obvious that the ATE leads to underestimation of mutations rates. In the first case (equal mutations rates), the bias in estimating *u _{i}* increases with

*i*and

*u*

_{4}is about two-thirds of the true value. MSEs of the estimates also suggest that the AII performs considerably better (except for

*u*

_{1}for which there is little difference between the two methods). For the second case, the downward bias in the estimates by the ATE is also obvious in all

*u*and the MSEs by ATE are appreciably larger than those by the AII. Another shortcoming of the ATE is that due to differential degrees of underestimation of

_{i}*u*, it can lead to rejection of certain hypotheses more often than specified by the given nominal level of significance. For example, for testing the hypothesis

_{i}*u*

_{1}=

*u*

_{4}, the ATE in the first case leads to nearly 12% rejection while the AII has <5% rejection at the 5% significance level. These results agree well with an earlier conclusion made from Table 2, which is that when the mean mutation rate is >10

^{−4}, the ATE starts to lose accuracy.

## Discussion

This article presents a significantly improved framework for statistical inference of germline mutation rates, with specific reference to *D. melanogaster*. This framework includes coalescent theory and an improved algorithm for simulating sample genealogies to obtain various coefficients, a method for computing the probabilities of mutation patterns, and a likelihood method for estimating mutation rates and testing hypotheses about the pattern of mutation rates. Statistical properties of the inference framework were investigated through simulation. The new approximation method for computing the probabilities of mutation patterns is more accurate than the previous method by Gao *et al.* (2011), particularly when mutation rates are high. Nevertheless, the previous method is sufficiently accurate for the data reported by Gao *et al.* (2011), and thus all major conclusions remain intact. The new likelihood-based inference exhibits desirable and expected properties, including reduced bias and smaller standard deviation with increasing number of families. The asymptotic χ^{2}-distribution for the likelihood-ratio test is sufficiently accurate when the number of families is reasonably large and for the sample size reported in Gao *et al.* (2011).

This theoretical study paves the way for analysis of data from families with three and four mutations. Furthermore, the theoretical framework reported here can be adapted for studying germline mutational distribution in other organisms and for analyzing data generated through DNA typing or sequencing sperm samples. To apply the framework to other organisms the nature and mode of cell propagation of their germline populations would need to be determined. Application to data generated by DNA typing or sequencing will likely have its own issues, such as data accuracy, but since it is more economical to sequence larger regions with fewer samples than shorter regions with larger samples, observing multiple mutations will likely be the norm. Therefore, the statistical framework of inference described here will be relevant.

## Acknowledgments

I thank Sara Barton for her editorial assistance. This work was partly supported by grants from the Chinese National Science Foundation (30570248 and 91231120 YF) and by the Betty Wheless Trotter Endowment Fund from The University of Texas Health Science Center.

## Appendix

### Asymptotic Covariance of the Maximum-Likelihood Estimates

The asymptotic covariance matrix of the maximum-likelihood estimates is the inverse of the following matrix:

Sinceit follows that (A1)(A2)Since *S* is of the formwhere *t*(*i*_{1}, *i*_{2}, *…* , *i _{k}*) is constant, it follows that (A3)(A4)Furthermore, let

*n*be the number of

_{i}*i*in

*i*

_{1}, … ,

*i*; then (A5) (A6)(A7)Putting the results of Equations A3–A7 into Equation A2, together with

_{k}**u**replaced by , will lead to the numerical value of .

## Footnotes

*Communicating editor: Y. S. Song*

- Received March 22, 2013.
- Accepted May 5, 2013.

- Copyright © 2013 by the Genetics Society of America