A Bayesian approach for correcting for partial plating in fluctuation experiments

Summary

The fluctuation experiment is the preferred method for estimating microbial mutation rates. A difficult task facing the data analyst is to infer the mean number of mutations from the number of mutant cells that only indirectly reflects the number of mutations. Partial plating, commonly practised in the laboratory, renders this task even more challenging by allowing only a portion of the mutant cells to be counted. In this paper, we propose a Bayesian approach to correcting for partial plating in the analysis of fluctuation experiments.

1. Introduction

Microbial mutation rates are routinely estimated in the laboratory using the fluctuation experiment devised by Luria & Delbrück (1943). A fluctuation experiment consists of a number of test tubes; each test tube contains a liquid culture, into which a small number (N ₀) of wild-type cells are seeded. In an incubation period, wild-type cells in the tubes freely divide and undergo the mutation of interest at an unknown rate. Each mutation engenders a mutant cell, and mutant cells freely divide – producing mutant daughter cells. At the end of the incubation period, the average number of wild-type cells per tube (N _T) is determined by titration. Meanwhile, the number of mutant cells existing in a tube at the end of the incubation period is determined by transferring the contents of that tube onto a solid, selective culture contained in a dish. This procedure is called plating. Through plating, wild-type cells are eliminated by a selective agent, but mutant cells proliferate and form visible colonies on the solid culture. The number of mutant cells in each tube prior to plating is determined by counting mutant colonies. A more detailed account of the fluctuation experiment is given in Foster (2006).

Of biological interest is the mean number of mutations that occur in a test tube during the incubation period, as this quantity is equivalent to the all-important mutation rate. Since each mutation gives rise to a first-generation mutant cell, which can spawn a large number of mutant offspring, in general, the number of mutant cells exceeds the number of mutations. The task facing the data analyst is to infer the mean number of mutations per test tube (customarily denoted by the symbol m) from the number of mutant cells, because, unlike mutant cells, mutations are not directly observable. This inference task is effected by means of so-called mutant distributions. A mutant distribution models the number of mutant cells in terms of the fundamental parameter m or some similar parameter. A well-known mutant distribution is the Luria–Delbrück distribution, obtained by Lea & Coulson (1949). This distribution, denoted by LD(m), is defined by the probability generating function (p.g.f.)

(1)

Ma et al. (1992) derived a recursive algorithm for computing the probability function of the distribution determined by this p.g.f. An algorithm for computing the maximum likelihood estimate (m.l.e.) of m using the Newton–Raphson procedure was derived by Zheng (2005), and is available in the software package SALVADOR 2.3 (Zheng, 2008b). Other more elementary methods for estimating m have been catalogued by Rosche & Foster (2000).

Inference about the parameter m is sometimes hindered by so-called imperfect plating efficiency, which refers to the phenomenon that a mutant cell existing in a tube succeeds in forming a visible colony with a probability ε<1. We call ε the plating efficiency. Imperfect plating efficiency can be caused by deliberately transferring only a proportion of a tube's contents to a dish during plating or by a cell's inability to survive the physical trauma of plating or by a combination of both. The former cause of imperfect plating efficiency is called partial plating. As microbial cells are highly efficient in forming colonies after being plated, it is justifiable to focus on partial plating when correcting for plating efficiency in estimating microbial mutation rates. Partial plating dates back to Luria & Delbrück (1943), who performed partial plating in 9 out of the 10 experiments reported in their classic paper. Today partial plating is still commonly practised because the requirement to plate the entire cultures can sometimes pose a formidable logistic challenge to the experimentalist. The necessity to relax this requirement by partial plating is well explained in Foster (2006, p. 200):

… this requirement restricts the volume of culture that can be used without concentrating the cells (which can be tedious with many cultures). However, if several mutant phenotypes are to be assayed in the same cultures, sampling is unavoidable. In addition, because a large culture contains more “information” than a small culture, it is better to plate a small aliquot from a large culture than all of a small culture if a proper correction can be applied …

The latest method for such a correction is a likelihood-based estimator of m proposed in Zheng (2008a). However, all exiting correction methods rely on the assumption that the precise value of ε is known. In practice, the actual plating efficiency is likely to be slightly different from the intended plating efficiency, however meticulous the experimentalist may be. To appreciate the consequence of this assumption, consider one of the two recent fluctuation experiments performed by Crane et al. (1996) to study a mutation of Salmonella Typhimurium. Each of these two experiments consists of 11 2·0-ml liquid cultures, but only a portion of 0·2 ml from each culture was plated for the detection of mutant colonies. We shall focus on the experiment that produced the following data.

(2)

In principle, from the above data maximum likelihood estimates of m and ε can be computed simultaneously. However, in practice, the joint likelihood function often causes numerical instability (Angerer, 2001). Even if one could overcome the numerical instability problem, the classic likelihood approach still does not seem to be the optimal way of estimating mutation rates in the setting of partial plating. In such a setting, the experimentalist intends to plate a prescribed portion ε₀ of contents from each tube. The joint likelihood function approach ignores this knowledge about ε₀ and is thus unacceptable from a philosophical perspective. One strategy for accounting for this knowledge is to maximize the likelihood function along the line ε=ε₀ in the mε-plane. Zheng (2008a) gives a Newton–Raphson type method to implement this strategy. In the experiment whose data are given in eqn (2), ε₀=0·1. Using SALVADOR 2.3, one finds that the m.l.e. of m is . However, we also note that the m.l.e. of m is extremely sensitive to changes in the presumed plating efficiency. Assuming the plating efficiency to be 0·08, one obtains an m.l.e. of 343·04. Changing the plating efficiency to 0·12, one has an m.l.e. of 243·42. Thus, the second strategy seems to be too restrictive. Here we propose a Bayesian middle ground between the two extreme strategies. We postulate that the actual plating efficiency may not be precisely ε₀, but is instead a realization of a random variable centred at ε₀ with a relatively small variance. Thus, two experiments with the same intended plating efficiency ε₀ may have different actual plating efficiencies due to different laboratory conditions, different experimentalists conducting the experiments or some other factors. Note that Asteris & Sarkar (1996) proposed a Bayesian approach to estimating m under the assumption of perfect plating efficiency.

2. Bayesian formulation and Markov chain simulation

A natural way to model imperfect plating efficiency is to view the plating procedure as a binomial thinning process applied to an LD(m) distribution (Armitage, 1952). Thus, the p.g.f. of the number of mutant colonies can be obtained by replacing z in eqn (1) with 1−ε+εz, which gives the p.g.f.

(3)

A distribution having the above p.g.f. is called an LD⁻(m, ε) distribution. (The latest work on this distribution is by Gerrish (2008).) If Y~LD⁻(m, ε), the probability mass function p(k; m, ε)=Prob(Y=k) does not have a closed-form expression. Using results obtained by Stewart (1991), Jones (1993a, b, 1994) and Angerer (2001), Zheng (2008a) gives the following recursive algorithm for computing p(k; m, ε) for k=0, 1, … , K. First, if ε⩽0·5, computing the following auxiliary sequence:

(4)

A variation of the above equations is recommended for ε>0·5 (Zheng, 2008a), but we need not consider such cases in the context of partial plating. Next, compute the probabilities recursively as follows:

(5)

The above algorithm permits relatively fast calculations of the joint likelihood function of m and ε.

Now let y ₁, … , y _n be the number of mutant colonies observed in an n-tube fluctuation experiment in which an ε-portion of liquid culture from each tube is plated. Thus, the y _i are an iid sample from an LD⁻(m, ε) distribution. The likelihood function is of the form

(6)

where p(·; m, ε) is the probability mass function of an LD⁻(m, ε) distribution. The likelihood function can be computed by first calculating p(k; m, ε) for k=0, 1, … , max_1⩽j⩽n y _j recursively using eqns (4) and (5) and then forming the product in eqn (6) using relevant terms. We relied on SALVADOR 2.3 (Zheng, 2008b) to compute the probability mass function.

An important issue in modelling the data is the choice of a prior distribution for the plating efficiency ε. As cell death caused by plating is assumed to be negligible, it appears a basic requirement that the prior concentrate in a relatively small neighborhood of ε₀ and be symmetrical about ε₀. The beta distribution is widely used to model proportions. One strategy is to match the mean; that is, let ε have a Beta(α, β) distribution with α/(α+β)=ε₀. Another strategy is to match the mode; that is, use a Beta distribution with (α−1)/(α+β−2)=ε₀. However, neither approach produces a distribution that is strictly symmetrical about ε₀. To meet the basic requirement, we depart from convention and choose a truncated normal distribution as a prior for ε. We employ the notation N _trunc(ε₀, υ, a, b) to denote a distribution whose density function has the form

We further let b−ε₀=ε₀−a to ensure symmetry of the truncated normal density function.

To complete the prior specification, we assign to the parameter m a relatively diffuse prior distribution of the form

(7)

The median of this distribution is m ₀ (Albert, 2007, p. 146). One can assign an empirical value to m ₀. For example, an accepted representative mutation rate for Escherichia coli cells is 10⁻⁹ per cell division. Therefore, for simplicity, one can set m ₀=10⁻⁹N _T. A sensible value for m ₀ can also be chosen by considering independent but similar existing experiments. From such an independent experiment one can derive several estimates of m, e.g., the m.l.e. of m, that can be used as parameter values for m ₀ in the prior density in eqn (7). A third option, which appears most useful in practice, is to derive a value for m ₀ directly from the data to be analysed. For example, one can easily obtain a candidate for m ₀ by utilizing Jones' formula

(8)

where is the median of the data (Jones et al., 1994). By replacing ε with the intended plating efficiency ε₀=0·1 in the above equation and applying the resulting formula to the data in eqn (2), we arrive at a sensible value for m ₀ which we denote by . This method of assigning a value to the parameter m ₀ in the prior in eqn (7) was inspired by the empirical Bayes idea (e.g. Carlin & Louis, 2009, p. 225). We found that posterior inference about m is not sensitive to the choice of a value for m ₀.

Following convention, we further assume that the prior distributions of ε and m are independent. Consequently, the posterior joint density of m and ε is proportional to

To facilitate posterior simulation, we make a change of variables

(9)

which transforms the parameter space onto the entire Euclidean plane. The Jacobian of this transform is

Accounting for the Jacobian of this transformation, the logarithm of the posterior density is, up to a constant term independent of θ₁ and θ₂,

(10)

where

Posterior simulation can be further facilitated by reducing the correlation between θ₁ and θ₂ via the following linear transform:

(11)

Here c is a real constant to be chosen by trial and error. Since the Jacobian of the transform in eqn (11) is a constant, the logarithm of the posterior density for ζ₁ and ζ₂ is, up to an additive constant, given by

(12)

Our inference about the fundamental parameter m is based on the log posterior joint density of ζ₁ and ζ₂ given in eqn (12). Using eqns (4) and (5), we adapted the standard ‘Metropolis within Gibbs’ method as outlined in the appendix to perform posterior inference. A simple way of choosing initial values for this algorithm is to set and ε=ε₀ in eqn (9). Since log((ε₀−a)/(b−ε₀))=0, it follows from eqn (11) that the initial states of the chain are

(13)

3. Illustration and discussion

We illustrate the above algorithm by analysing the data given in eqn (2). We first choose a truncated normal distribution as a prior for m. Recall that the experimentalist intended to plate 10% of the contents from each tube, which dictates that we set ε₀=0·1. We assume that the actual plating efficiency falls within 50% of either side of the intended plating efficiency, and hence we set a=0·05 and b=0·15. The choice of a value for the variance parameter υ can reflect the experimentalist's uncertainty about the actual plating efficiency, which can be guided by a simple graphical approach. By inspection of plots of various truncated normal density functions like Fig. 1, we choose υ=10⁻⁴ to reflect our uncertainty about the actual plating efficiency. To get a feel for this particular prior, N _trunc(0·1, 10⁻⁴, 0·05, 0·15), we calculated the following two prior probabilities: P(0·09<ε<0·11)= 0·68 and P(0·08<ε<0·12)=0·95. The second probability implies that the experimentalist is 95% certain that the actual plating efficiency falls within 20% of either side of the intended plating efficiency. We next focused attention on the choice of the linear transform constant c appearing in eqn (11) and the two scale parameters, c ₁ and c ₂, defined in the appendix. For each trial triplet (c, c ₁, c ₂), we ran the algorithm for a few thousands of iterations to examine the two acceptance rates and the sample correlation coefficient, ρ=corr(ζ₁, ζ₂). We aimed at obtaining acceptance rates of around 0·44 (see, e.g. Albert, 2007, p. 105) and a correlation ρ of relatively small magnitude. In addition, we examined contour plots, trace plots and autocorrelation plots to assess convergence of the chain. We finally chose c=2·15, c ₁=0·30 and c ₂=0·75. We then started our ‘final’ Markov chain by setting in eqn (10) and using eqn (13) with the same value. We ran the algorithm outlined in the appendix for 50 000 iterations after a burn-in of 5000 iterations. Every 10th iterate is saved for summaries and inferences (see Figs. 2–4 for more details about these saved draws). This recipe yielded ρ=−0·013 and acceptance rates of 0·44 and 0·43 for ζ₁ and ζ₂, respectively. Inference about m is based on the induced chain of e^ζ₁. The mean of this induced chain is and the median is 282·3. Both are close to the m.l.e., , obtained using the existing method of Zheng (2008a). Since the present approach takes into account a degree of uncertainty about the actual plating efficiency, the 95% credible interval for m, which is (218·2, 359·7), is wider than the corresponding confidence interval (Zheng, 2008b), which is (232·5, 332·4).

Fig. 1. Three truncated normal density functions N _trunc(0·1, υ, 0·05, 0·15) with different values of υ, which can be considered as prior distributions for ε when the intended plating efficiency is ε₀=0·1.

Fig. 2. Contour plot of the posterior joint density of ζ₁ and ζ₂ along with the last 500 draws of a simulated chain for the data given in eqn (2).

Fig. 3. Trace plot of the last 2000 ζ₁ draws for the analysis of data given in eqn (2).

Fig. 4. Histogram of 5000 simulated values of the mean number of mutations m=e^ζ₁. In this simulation for the data in eqn (2) the prior for ε was N _trunc(0·1, 10⁻⁴, 0·05, 0·15).

Although the proposed method cannot be automated to the extent to which the maximum likelihood method can, the method is still simple enough for routine practical use. Sensible choice of the three tuning parameters, c, c ₁ and c ₂, might require time and experience on the part of the data analyst, but it is worth noting that the above simulation recipe, i.e., the length of a burn-in period and the number of iterations after the burn-in period, can tolerate a relatively wide range of values of the tuning parameters. As an illustration, we re-ran the above simulation for 1000 times. Each time we followed the same recipe, but allowed the three tuning parameters to be determined by the following equations:

where r ₁, r ₂ and r ₃ were drawn from Uniform(0,1). The 0·025, 0·500 and 0·975 posterior quantiles are summarized in Fig. 5. The box plots clearly indicate that all the simulations led essentially to the same inference about m.

Fig. 5. The same simulation recipe described in the last section were repeated 1000 times with random tuning parameters. From the top panel to the bottom, the box plots describe the 0·025, 0·500 and 0·975 posterior quantiles, respectively, for the data given in eqn (2).

The proposed approach makes explicit a subjective ingredient in the estimator of m, which was hidden in all non-Bayesian estimators as the assumption of absolute certainty about actual plating efficiency is also subjective. The new approach offers a practical recipe for obtaining a better indication of the range of m that can be expected in view of the data at hand.