The Theoretical Basis of qPCR and ddPCR Copy Number Estimates: A Critical Review and Exposition

1. Introduction

The polymerase chain reaction (PCR) is a fundamental tool of molecular biology that has found widespread application in the biological and environmental sciences. An early form of PCR was developed in the late 1960s by Kjell Kleppe, who called it “repair replication” [1]. The modern form was conceived in 1983 by Kary Mullis, who subsequently developed a practical and efficient method that he published in a series of papers with various collaborators [2,3,4]. Mullis shared the 1993 Nobel Prize in Chemistry for this work.

PCR is essentially a molecular copy machine that can rapidly generate billions of copies of a short DNA sequence. In environmental studies, it can be used in conjunction with gel electrophoresis and a fluorescent DNA stain to provide a visual presence/absence assay for target DNA sequences from specific groups or organisms, e.g., [5,6]. Standard PCR, however, provides no quantitative information about the concentrations of the target sequences in a sample. For example, while standard PCR can detect target sequences characteristic of enterococci fecal fecal indicator bacteria in recreational beach water samples, it cannot determine their concentration and therefore cannot be used to either establish or assess compliance with a numerical criterion for protecting human health, such as a Beach Action Value [7]. Various elaborations of PCR have been developed since the 1980s to address this important limitation. Of these quantitative PCR methods, the two that are most widely used at the present time are real-time quantitative PCR (qPCR) and droplet digital PCR (ddPCR).

Real-time quantitative PCR was first developed by Higuchi and others in the early 1990s [8]. The key idea underlying its ability to determine the concentration of a target sequence in a sample is to introduce a unique fluorescent label or dye whose level of fluorescence is tightly linked to amplification of the target DNA sequence by PCR, so that fluorescence of a sample increases as the number of amplicons increases. Unlike the concentration of the amplicons themselves, fluorescence can easily be monitored in real time and, with additional information from a standard curve, the temporal pattern of increase can be used to back-calculate the initial number of copies in the sample. The dependence of this approach on monitoring the polymerase chain reaction in real time is the key feature that distinguishes it from endpoint types of quantitative PCR, such as chamber-based digital PCR (cdPCR) and ddPCR [9]. The acronym “qPCR” is now generally recognized to denote specifically real-time quantitative PCR.

The original forms of digital PCR also were developed in the 1990s [10,11]. These methods were based on the idea of diluting extracted DNA containing the target sequence to such a degree that most wells in a multi-well PCR plate will contain either no copies of the sequence or only a small number. Sample wells containing the extracted DNA and an assay mix that included a fluorescent dye were processed simultaneously as replicate endpoint PCR reactors. The concentration of the target DNA sequence was estimated from the proportion of wells that failed to fluoresce above a user-specified threshold level (hence, contained no copies of the target sequence), using a probability argument based on the Poisson distribution. Important advantages over qPCR were reduced sensitivity to PCR inhibition and the fact that no standard curve was required.

The modern droplet-based form of digital PCR was developed during the late 1990s and early 2000s and employs microfluidics technology [12,13,14,15]. With this method, extracted DNA containing the target sequence (along with an assay mix that includes a fluorescent dye) is partitioned into thousands of microscopic droplets per sample well that are so small that most will contain either no copies of the target sequence or only a very small number. The temperature cycling that drives PCR amplification is then performed simultaneously on all the droplets, which serve as replicate endpoint PCR reactors. Each droplet fluoresces above a user-specified threshold level at the reaction endpoint if and only if at least one copy of the target sequence was included in the droplet when formed (and was then amplified). As in the dilution form of digital PCR, the proportion of droplets that fail to fluoresce above the threshold, as determined by a droplet reader, is used to estimate the concentration of the target sequence, and no standard curve is required. Thus, ddPCR achieves the same outcome as earlier forms of digital PCR but with far less dilution of the sample and a vastly larger number of replicate reactions. ddPCR currently is more expensive than qPCR and has a somewhat longer sample processing time, but it also has important advantages in environmental applications (notably, no dependence on a standard curve and reduced sensitivity to interference) and therefore is increasingly being used in place of qPCR.

The purpose of this paper is to systematically develop what we believe to be the mathematical and statistical foundations of qPCR and ddPCR in sufficient detail to identify key underlying assumptions that may affect the accuracy and precision of concentration estimates and, more broadly, to answer two fundamental questions to the reader’s and our own satisfaction:

• How does qPCR make it possible to estimate the initial copy number in an environmental sample by monitoring, in real time, the increasing fluorescence during successive PCR cycles, using additional information derived from a standard curve?
• How does ddPCR make it possible to estimate the initial copy number in a sample by determining the proportion of droplets that do not fluoresce above background at the reaction endpoint, without requiring a standard curve?

Some of the details of our presentation will no doubt reflect the specific analytical instruments we normally use for qPCR and ddPCR analysis, which are the Applied ${\mathrm{Biosystems}}^{\circledR }{\mathrm{StepOnePlus}}^{\mathrm{TM}}$ real-time qPCR system and the ${\mathrm{Bio-Rad}}^{\circledR }{\mathrm{QX200}}^{\mathrm{TM}}$ ddPCR system. The basic principles, however, appear to be general.

Software accompanying laboratory instruments that perform qPCR and ddPCR carry out all the necessary calculations and report the estimated copy numbers, but we are philosophically opposed to trusting proprietary software to perform calculations that are not fully documented and whose software implementation cannot be checked. Some of the necessary information about these calculations is available in the literature, e.g., [9,16,17,18] and on instrument manufacturer and vendor websites, e.g., [19,20], but all the accounts we are aware of are incomplete, and some of the information is overly simplified and not entirely correct. We therefore set out to determine for ourselves whether rigorous quantitative theories can be developed for qPCR and ddPCR that justify placing trust in the numbers that commercial software produces. As will be shown, the answer is a qualified “yes” for both qPCR and ddPCR.

The paper is organized as follows. We address qPCR first, then ddPCR. For each method, we begin with an overview, then develop the underlying mathematical and statistical theory, and finally present an example where we apply the theory to some of our own data. We conclude with a brief discussion comparing what we see as the strengths and weaknesses of the two methods when used to quantify DNA target sequences in environmental samples, especially in studies of aquatic systems.

2. Real-time Quantitative PCR

2.1. A Brief Overview of qPCR

Numerous protocols for qPCR analysis of environmental samples are available in the literature and on instrument manufacturer and vendor websites (an example of a qPCR workflow that has been widely used in the state of Michigan for monitoring E. coli contamination at recreational beaches is presented in Appendix A). Regardless of which protocol one uses, the first step in this process is sample preparation. In studies quantifying the abundances of aquatic microorganisms such as bacteria or unicellular algae, cells of the focal taxonomic group (as well as others) are collected from water samples, usually by filtration, and then physically or chemically lysed, releasing the genetic material they contain for subsequent exposure to the reagents employed in qPCR analysis.

Crude extracts from physical lysis procedures such as bead milling, consisting of the supernatants of these samples after centrifugation, are often directly analyzed with qPCR. However, this procedure is not very effective at removing common PCR inhibitors such as tannins, chlorophyll, and humic/fulvic acids that are collected with the cells, and it is not applicable with most chemical lysis procedures. Therefore, sample preparation processes will often include a DNA purification procedure to remove contaminants and, if necessary, concentrate the recovered genetic material.

A variety of procedures can be used to purify and concentrate extracted samples. The most common types involve passage of the crude lysates through silica membranes, where a salt gradient binds all nucleic acids in the sample to the membrane [21,22]. After the membrane is washed to remove inhibitors and other contaminants, the genetic material is eluted and is ready for amplification.

Like standard (nonquantitative) PCR, qPCR creates a large number of copies of (i.e., amplifies) a target DNA sequence in a prepared sample if that sequence is present. Unlike standard PCR, however, it also quantifies the number of copies of the target sequence that were present before amplification. Quantification is achieved by combining the prepared sample with an assay mix (which contains forward and reverse primers that bind to the endpoints of the target DNA sequence, a fluorescent DNA dye or probe that may or may not be specific to the target sequence, a reference dye (Section 2.3.1), deoxynucleotide triphosphates (dNTPs), DNA polymerase, salts, and a buffer), then measuring in real time the increasing intensity of fluorescence that occurs during a series of consecutive amplification cycles. Two alternative types of fluorescent dye can be used, with their most common representatives being ${\mathrm{TaqMan}}^{\circledR }$ probes and the ${\mathrm{SYBR}}^{\circledR }$ Green DNA stain. TaqMan probes are homologous to a segment of the DNA target sequence, contain separate fluorescent and fluorescence-quenching moieties (the reporter and quencher, respectively), and exhibit little fluorescence until the reporter is cleaved from the probe by the Taq DNA polymerase used to make copies of the target sequence, physically separating it from the quencher. By contrast, SYBR Green is a nonspecific DNA stain that binds mainly to the minor groove of any double-stranded DNA that is present [23] and fluoresces little unless bound to double-stranded DNA. To condense our exposition, we will restrict attention throughout this paper to TaqMan probes, which are more specific than SYBR Green and can be used for multiplexing.

Each PCR amplification cycle consists of three phases that are initiated and terminated by temperature changes, with the specific temperatures depending on the template DNA: (1) DNA denaturation to separate the two strands of DNA (temperature: typically 94–98° C), (2) annealing of sequence-specific forward and reverse primers to the separate strands, and of the sequence-specific TaqMan probe to the appropriate single strand (temperature: typically 50–70°C), and (3) extension of both primers by a thermostable DNA polymerase (temperature: typically 68–72°C). During the extension phase, DNA polymerase cleaves the fluorescent moiety of the probe attached to the target sequence from the quencher moiety, producing a single fluorescing free reporter for each copy of the target DNA sequence that is present. As the cycles repeat, the number of copies of the target sequence increases geometrically (Section 2.2.1) while the number of fluorescing free reporters increases approximately geometrically (Section 2.2.2). The corresponding fluorescence level due to free reporters is initially much lower than the relatively constant background level and is therefore masked by background fluorescence (Figure 1). After several cycles (the actual number corresponds closely to the initial concentration of the target sequence), the rapidly increasing fluorescence level due to free reporters exceeds the background level and becomes measurable, with log-transformed fluorescence exhibiting an approximately linear increase as cycle number increases (Figure 1, Figure 2). Eventually the residual concentration of one or more PCR reagents is reduced to levels that begin to noticeably depress the rate of increase in fluorescence [24]. Fluorescence then exhibits a nonlinear decelerating pattern of increase on a logarithmic scale and approaches an asymptote that generally corresponds to the endpoint of standard PCR.

For reasons that will become clear in Section 2.2, estimates of the number or concentration of target sequence copies (TSC; we will use “TSC” as singular or plural, depending on context) in the original sample must be obtained from the linear portion of the log-fluorescence curve. Intuitively, this is the range of amplification cycles where the geometric increase in copy number, as reflected in the approximately geometric increase in fluorescence, reveals itself: the fluorescence “signal” due to amplification has increased sufficiently to separate itself from background “noise” but depletion of PCR reactants is not yet sufficient to cause a clear deviation from geometric increase. A threshold level of fluorescence whose log-transformed value lies within this linear region is chosen, and the fractional cycle number at which fluorescence in a given sample crosses the threshold (i.e., the threshold cycle) is estimated by interpolation. The fundamental idea on which qPCR quantification rests is very simple: the smaller the number of amplification cycles required for fluorescence to cross this threshold, the more copies of the target sequence must have been present in the original sample.

2.2. Basic Equations of the qPCR Process

2.2.1. Number of Target Sequence Copies

The schematic diagram in Figure 3 illustrates the increase in the number of copies of a target DNA sequence and the corresponding increase in the number of fluorescing free reporters that occur as the number of PCR cycles increases. Each cycle consists of three main temperature-dependent steps (the specific temperatures listed here and in Figure 3 are commonly used for E. coli DNA):

• Denaturation of DNA at 95°C
• Annealing of target-specific forward and reverse primers and probes to the target sequence at 50°C
• Extension of primers by DNA polymerase and concomitant cleavage of probes at 72°C to yield free reporters that fluoresce when excited with the appropriate wavelength of light.

For convenience, we will refer to this cycle as the DAE (Denaturation, Annealing, Extension) cycle.

Figure 3 also shows that if every target sequence is amplified during each DAE cycle, then each TSC present at the beginning of a given cycle produces one new TSC and one new free reporter during the cycle. TSC and free reporters already present at the beginning of a cycle persist, so the total number of each at the end of the cycle is the sum of the old and new entities. In practice, some of the TSC present at the beginning of a cycle fail to produce a new copy during the cycle, so the process is not 100% efficient. We will assume, however, that each new TSC produced is accompanied by a new free reporter.

Assuming no degradation of existing TSC occurs, the number of copies present at the end of c DAE cycles is related to the number present at the end of $c-1$ cycles by the numerical balance law,

$$\left[\begin{array}{c}\mathrm{Number}\mathrm{of}\\ \phantom{,}\mathrm{TSC}\mathrm{after}\phantom{,}\\ c\mathrm{cycles}\end{array}\right]=\left[\begin{array}{c}\mathrm{Number}\mathrm{of}\\ \phantom{,}\mathrm{TSC}\mathrm{after}\phantom{,}\\ c-1\mathrm{cycles}\end{array}\right]+\left[\begin{array}{c}\mathrm{Number}\mathrm{of}\mathrm{new}\\ \phantom{,}\mathrm{TSC}\mathrm{produced}\phantom{,}\\ \mathrm{during}\mathrm{cycle}c\end{array}\right],c=1,2,3,\dots$$

(1)

To translate this verbal statement into mathematical symbols, let $x\left(c\right)$ be a dimensionless function mapping cycle number c to the number of TSC present at the end of the cycle, and let $\epsilon \in (0,1]$ denote the average proportion of TSC that are amplified during any given DAE cycle. We follow EPA Draft Method C [25,26] in expressing copy numbers as numbers per reaction (which is dimensionless) instead of numbers per volume (which has dimension ${\mathrm{Length}}^{-3}$), with the reaction volume unspecified but assumed consistent in all equations. Then $\epsilon$ represents the proportional efficiency of amplification. If amplification is fully efficient ($\epsilon =1$), each TSC present at the beginning of a DAE cycle produces one new copy during the cycle. In the more realistic case where amplification is less than fully efficient ($0<\epsilon <1$), most existing copies will produce new copies but some will not. In this case, each TSC present at the beginning of DAE cycle c will produce an average of $\epsilon$ new copies during the cycle, so the total number of new TSC produced during the cycle will be given by

$$\left[\begin{array}{c}\mathrm{Number}\mathrm{of}\mathrm{new}\\ \phantom{,}\mathrm{TSC}\mathrm{produced}\phantom{,}\\ \mathrm{during}\mathrm{cycle}c\end{array}\right]=\epsilon x(c-1).$$

Replacing the third term of the balance law in Eq. (1) with this expression, and replacing the first and second terms with $x\left(c\right)$ and $x(c-1)$, we find that

$$\begin{array}{ccc}\hfill x\left(c\right)& =\hfill & x(c-1)+\epsilon x(c-1)\hfill \\ & =\hfill & (1+\epsilon )x(c-1)\hfill \\ & =\hfill & \lambda x(c-1),c=1,2,3,\dots ,\hfill \end{array}$$

(2)

where $\lambda$ is the amplification factor defined by

$$\lambda :=1+\epsilon >1.$$

(3)

Symbols and associated dimensions used in the qPCR equations here and below are summarized in Table 1.

Eq. (2) is a homogeneous first-order linear difference equation with initial condition $x\left(0\right)=x_0>0$, where initial copy number $x_0$ is treated as a constant. Solving this equation by backward iteration, we find that

$$x\left(c\right)=x_0{\lambda }^c,c=0,1,2,\dots ,$$

(4)

which is geometric in cycle number c (cf. [18]). (Note: $x\left(c\right)$ is a geometric sequence, because $x(c+1)/x\left(c\right)=\mathrm{constant}=\lambda$ for all c; it is also exponential in sequence index c.) Taking base-10 logarithms on both sides (any base will do, but base 10 is customary in the qPCR literature) then yields

$${log}_{10}\left(x\left(c\right)\right)={log}_{10}\left(x_0\right)+{log}_{10}\left(\lambda \right)c.$$

(5)

Thus, the logarithm of the number of copies present after c DAE cycles is a linear function of the number of cycles, with intercept ${log}_{10}\left(x_0\right)$ and slope ${log}_{10}\left(\lambda \right)>0$. The slope is mainly important in assessing efficiency of the qPCR process, while the intercept is important in estimating the initial copy number (see below).

2.2.2. Number of Free Reporters

As noted in Section 2.1, measuring the fluorescence of free reporters is how qPCR quantifies the number of TSC, exploiting the intimate relationship between the production of new TSC and cleavage of annealed TaqMan probes during primer extension. It is therefore necessary to develop an equation for the number of free reporters as a function of DAE cycle c, and then to connect that equation to fluorescence.

Let $r\left(c\right)$ denote the number of free reporters at the end of DAE cycle c. Each TSC present at the beginning of cycle c (hence, each TSC present at the end of cycle $c-1$) produces, on average, $\epsilon$ new reporters (one for each new TSC) by the end of the cycle. Therefore,

$$\begin{array}{ccc}\hfill r\left(c\right)& =\hfill & r(c-1)+\epsilon x(c-1)\\ & =\hfill & r(c-1)+\epsilon x_0{\lambda }^{c-1}.\end{array}$$

(6)

This recurrence relation holds for each pair of consecutive DAE cycles. Iterating backward, and requiring $r\left(0\right)=0$ (because no DAE cycles have occurred yet when $c=0$), we find that

$$\begin{array}{}& \hfill r\left(c\right)& =\hfill & r\left(0\right)+\epsilon x_0(1+\lambda +{\lambda }^2+\cdots +{\lambda }^{c-1})& & =\hfill & \epsilon x_0({\lambda }^c-1)/(\lambda -1)=\epsilon x_0({\lambda }^c-1)/\epsilon \mathrm{(7)}& & =\hfill & x_0({\lambda }^c-1)\mathrm{(8)}& & =\hfill & x\left(c\right)-x_0,c=0,1,2,\dots \end{array}$$

Noting that $x\left(c\right)-x_0$ is the number of new TSC created by amplification during DAE cycles 1 through c, the last of these equations states that the cumulative number of free reporters present after c cycles is equal to the cumulative number of new TSC that have been created. Intuitively, this must be the case, because the primer extension that produces each new TSC also produces one new free reporter.

2.2.3. Fluorescence

Assume that the fluorescence intensity $f\left(c\right)$ of a sample of fixed volume at the end of DAE cycle c is proportional to the number $r\left(c\right)$ of free reporters and hence, for a fixed reaction volume, to the concentration of free reporters. Then

$$f\left(c\right)=\phi r\left(c\right),$$

(9)

where $\phi$ is a multiplicative constant that converts the number of free reporters (in a fixed volume) to fluorescence intensity. It follows from Eq. (7) that

$$f\left(c\right)=\phi x_0\left({\lambda }^c-1\right).$$

(10)

This equation gives the fluorescence intensity at the end of a cycle as a function of cycle number c. It will allow us to estimate the initial copy number $x_0$, once we have adapted it to address certain complications that arise in the laboratory.

2.3. Adapting the Basic Equations for Use in the Laboratory

When applying the theory outlined in the previous section to real samples, several complications arise. These include background fluorescence, well effects, and inhibition. The basic theory can be adapted to address all of these complications. Continual changes in laboratory instruments are likely to alter some of the details we discuss regarding how the basic qPCR equations are adapted to the laboratory, but it should always be possible to modify the mathematical and statistical theory we outline to accommodate these new details and thereby retain a sound theoretical basis for qPCR copy number and concentration estimates.

2.3.1. Accounting for Background Fluorescence and Well Effects

Background fluorescence at the emission wavelength of the reporter occurs in all plate wells containing samples (or standards) and master mix, mainly due to incomplete quenching of reporter fluorescence in the intact probe (e.g., [19,27]). Total reporter fluorescence in each well is the sum of background fluorescence and fluorescence of the free reporter that occurs after cleavage of the probe by DNA polymerase and exposure to light of the excitation wavelength. To express this idea in symbols, let $B_w$ denote the level of background fluorescence in well w, and suppose this background level remains approximately constant as the number of amplification cycles increases. Then the total intensity $g_w\left(c\right)$ of reporter fluorescence in well w is given by

$$g_w\left(c\right)=f\left(c\right)+B_w=\phi x_0({\lambda }^c-1)+B_w,$$

(11)

where we used Eq. (6).

Another complication that must be accounted for is that identical concentrations of a sample or standard in different plate wells typically show different measured fluorescence intensities at the emission wavelength of the reporter. This well effect has several potential causes, including bubbles, condensation, evaporation, and in some instruments, differences in the length of the light path (e.g., [27,28]). To account for a well effect, we assume that the measured intensity of total reporter fluorescence $R_w\left(c\right)$ in well w in cycle c is given by the product of the true level $g_w\left(c\right)$ and a dimensionless well-effect factor ${\kappa }_w>0$. The measured intensity of total reporter fluorescence in well w in cycle c is therefore given by

$$R_w\left(c\right)={\kappa }_w{g}_w\left(c\right)={\kappa }_w\left[\phi x_0({\lambda }^c-1)+B_w\right].$$

(12)

The combination of background fluorescence and well effect in measurements of fluorescence intensity can be removed by including a passive reference dye in the assay master mix that fluoresces at a different emission wavelength than the reporter dye (the most common reference dye is carboxy-X-rhodamine, or ROX). Let $\rho$ denote the fluorescence intensity of the passive reference dye in each well of the qPCR plate that contains the assay master mix, with $\rho$ assumed to be the same for all samples (or standards). The well effect modifies this fluorescence to the same extent that it modifies fluorescence of the reporter dye, because the cause is the same. Therefore, the measured fluorescence intensity $P_w$ of the reference dye in well w is given by

$$P_w={\kappa }_w\rho$$

(13)

for all cycles c. Dividing the measured fluorescence intensity of the reporter in cycle c, as given by Eq. (12), by the measured fluorescence level of the reference dye, as given by Eq. (13), we obtain the normalized reporter fluorescence $R{n}_w\left(c\right)$ in well w,

$$R{n}_w\left(c\right)={\displaystyle \frac{R_w\left(c\right)}{P_w}}={\displaystyle \frac{\phi }{\rho }}{x}_0\left({\lambda }^c-1\right)+{\displaystyle \frac{B_w}{\rho }}.$$

(14)

$R{n}_w\left(c\right)$ is dimensionless, because it is the ratio of two measures of fluorescence intensity (however, software for qPCR instruments may report $R{n}_w\left(c\right)$ in “relative fluorescence units”, or RFU).

In practice, fluorescence due to free reporters is so low during the first 15 or so DAE cycles that it would not be detectable even in the absence of background fluorescence. Therefore, the measured normalized fluorescence during these preliminary cycles is due entirely to background fluorescence, and the ratio $B_w/\rho$ can therefore be estimated (this is often done by averaging the normalized fluorescence measurements in each well over cycles 3 to 15). Assuming this estimate is in fact equal to $B_w/\rho$ and subtracting it from both sides of Eq. (14), we obtain the following equation for the difference, $\Delta Rn\left(c\right)$:

$$\Delta Rn\left(c\right)=R{n}_w\left(c\right)-{\displaystyle \frac{B_w}{\rho }}={\displaystyle \frac{\phi }{\rho }}{x}_0\left({\lambda }^c-1\right).$$

(15)

Note that $\Delta Rn\left(c\right)$ depends on initial copy number $x_0$ but that the well effect has been accounted for and background fluorescence has been removed. We may refer to $R{n}_w\left(c\right)$ as the normalized reporter fluorescence, to $B_w/\rho$ as the normalized background fluorescence, and to $\Delta Rn\left(c\right)$ as the background-corrected normalized reporter fluorescence. All three quantities are dimensionless.

2.3.2. An Equation for the Threshold Cycle

The threshold cycle $c_{\mathrm{T}}$ for a sample or standard is the fractional cycle at which background-corrected normalized reporter fluorescence $\Delta Rn\left(c\right)$ crosses a user-specified threshold level $\tau$ (some authors call the threshold cycle the “quantification cycle” and denote it by $C_q$, but this term disguises the direct meaning of the threshold cycle, which is simply the fractional cycle at which $\Delta Rn\left(c\right)$ crosses the threshold). In practice, $c_{\mathrm{T}}$ values typically range from roughly 20 to 35 cycles, and $\lambda$ is required to be very close to 2 (so $\epsilon$ is very close to 1) to ensure data quality.

By definition of the threshold cycle, we must have $\Delta Rn\left(c\right)=\tau$ when $c=c_{\mathrm{T}}$ in Eq. (15), temporarily relaxing the restriction of c to the integers. Thus, we require

$${\displaystyle \frac{\phi }{\rho }}{x}_0\left({\lambda }^{c_{\mathrm{T}}}-1\right)=\tau .$$

(16)

This equation can be rearranged to obtain

$${\lambda }^{c_{\mathrm{T}}}={\displaystyle \frac{\tau \rho }{\phi x_0}}+1={\displaystyle \frac{\tau \rho }{\phi x_0}}\left(1+{\displaystyle \frac{\phi x_0}{\tau \rho }}\right),$$

(17)

where both sides of the equation are dimensionless. Taking base-10 logarithms,

$$c_{\mathrm{T}}{log}_{10}\left(\lambda \right)={log}_{10}\left({\displaystyle \frac{\tau \rho }{\phi x_0}}\right)+{log}_{10}\left(1+{\displaystyle \frac{\phi x_0}{\tau \rho }}\right).$$

(18)

The last term on the right is negligible in practice and can be ignored (see Appendix B). To a very good approximation, then,

$$c_{\mathrm{T}}{log}_{10}\left(\lambda \right)={log}_{10}\left({\displaystyle \frac{\tau \rho }{\phi x_0}}\right).$$

(19)

Dividing both sides of this equation by ${log}_{10}\left(\lambda \right)>0$, we obtain an equation for $c_{\mathrm{T}}$ as a nonlinear function of initial concentration $x_0$:

$$\begin{array}{}\mathrm{(20)}& \hfill c_{\mathrm{T}}& =\hfill & {\displaystyle \frac{1}{{log}_{10}\left(\lambda \right)}}{log}_{10}\left({\displaystyle \frac{\tau \rho }{\phi x_0}}\right)& & =\hfill & {\displaystyle \frac{{log}_{10}(\tau \rho /\phi )}{{log}_{10}\left(\lambda \right)}}-{\displaystyle \frac{1}{{log}_{10}\left(\lambda \right)}}·{log}_{10}\left(x_0\right)\mathrm{(21)}& & =\hfill & {\beta }_0+{\beta }_1{log}_{10}\left(x_0\right)\end{array}$$

[cf. 18, p. 233], where $x_0\ge 1$, ${log}_{10}\left(x_0\right)\ge 0$, and

$$\begin{array}{ccc}\hfill {\beta }_0& =\hfill & {\displaystyle \frac{{log}_{10}(\tau \rho /\phi )}{{log}_{10}\left(\lambda \right)}}\end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\beta }_1& =\hfill & -{\displaystyle \frac{1}{{log}_{10}\left(\lambda \right)}}.\end{array}$$

(23)

Note that both sides of Eq. (21) are dimensionless, as are coefficients ${\beta }_0$ and ${\beta }_1$.

Eq. (21) can be fitted to a set of calibration data (measured $c_{\mathrm{T}}$ versus known ${log}_{10}\left(x_0\right)$ values for several different standards) to estimate coefficients ${\beta }_0$ and ${\beta }_1$, which in turn can be used to estimate amplification factor $\lambda$ and parameter ratio $\rho /\phi$ and their 95% confidence intervals. The statistical problem of how best to fit Eq. (21) to calibration data is nontrivial and is discussed below in Section 2.4. Recalling that $\lambda \approx 2$ in practice, it follows that ${log}_{10}\left(\lambda \right)>0$ and hence ${\beta }_1<0$, using Eq. (23). Thus, the threshold cycle is a linearly decreasing function of ${log}_{10}\left(x_0\right)$, with intercept ${\beta }_0>0$ and slope ${\beta }_1<0$ defined by Eqs. (22) and (23). This result implies that $c_{\mathrm{T}}={\beta }_0$ for $x_0=1$ and $c_{\mathrm{T}}<{\beta }_0$ for all $x_0>1$.

The fitted equation for $c_{\mathrm{T}}$ as a linear function of ${log}_{10}\left(x_0\right)$ is called the standard curve or calibration curve. As we discuss in Section 2.3.3, the standard curve can be rearranged to obtain an equation for ${log}_{10}\left(x_0\right)$ as a linear function of $c_{\mathrm{T}}$, which can then be used to estimate unknown concentrations in samples from measured $c_{\mathrm{T}}$ values. This procedure is called inverse prediction or classical calibration in the statistical literature.

An empirical example of the relationship between $c_{\mathrm{T}}$ and ${log}_{10}\left(x_0\right)$ is shown in Figure 4 (left panel). The data come from a qPCR calibration run employing five standards, with three replicates per standard. The dashed line is a fitted linear (ordinary) least-squares regression model. Note that the trend is linear to a reasonably good approximation.

2.3.3. Estimating the Initial Copy Number

Once parameters ${\beta }_0$ and ${\beta }_1$ have been estimated by fitting Eq. (21) to a set of calibration data (see Section 2.4), the inverted equation can be used to estimate the unknown initial number of TSC per reaction in field samples, based on their measured values of $c_{\mathrm{T}}$. Using the classical calibration approach to inverse prediction, we rearrange Eq. (21) and solve for ${log}_{10}\left(x_0\right)$ to obtain

$$\begin{array}{}\mathrm{(24)}& \hfill {log}_{10}\left(x_0\right)& =\hfill & (c_{\mathrm{T}}-{\beta }_0)/{\beta }_1\mathrm{(25)}& & =\hfill & ({\beta }_0-c_{\mathrm{T}})/\left|{\beta }_1\right|\mathrm{(26)}& & =\hfill & {log}_{10}(\tau \rho /\phi )-{log}_{10}\left(\lambda \right)c_{\mathrm{T}},\end{array}$$

which implies that

$$x_0={10}^{({\beta }_0-c_{\mathrm{T}})/\left|{\beta }_1\right|}={10}^{{log}_{10}(\tau \rho /\phi )-{log}_{10}\left(\lambda \right)c_{\mathrm{T}}}={\displaystyle \frac{\tau \rho }{\phi }}{\lambda }^{-c_{\mathrm{T}}}$$

(27)

[e.g., 29, p. 36]. (The last form on the right can be derived directly from Eq. (16) by noting that $1/({\lambda }^{c_{\mathrm{T}}}-1)$≈ ${\lambda }^{-c_{\mathrm{T}}}$, given that $\lambda \approx 2$ and $c_{\mathrm{T}}\ge 20$.)

Eq. (27) gives a point estimate of the initial copy number but gives no sense of its uncertainty. This uncertainty has two main components: one due to measurement error in estimates of the number of TSC per reaction in the standards used to fit Eq. (21), and the other due to statistical error in estimates of ${\beta }_0$ and ${\beta }_1$, given the standard estimates that were used.

To address the first source of error, suppose the set of standards are prepared by first preparing the highest standard and then preparing the others by a series of dilutions. Asssuming the dilutions are done properly, the main source of error will be in the estimated number of TSC per reaction in the highest standard. Let ${\tilde{x}}_{\max }$ and $x_{\max }$ denote the estimated and true numbers of TSC per reaction in the highest standard, and suppose the estimate contains proportional error $\eta x_{\max }$, so that

$${\tilde{x}}_{\max }=(1+\eta )x_{\max }.$$

(28)

The other standards will be estimated based on ${\tilde{x}}_{\max }$ and their dilution factors, so the estimate for the i-th standard will be

$${\tilde{x}}_i=D_i{\tilde{x}}_{\max }=(1+\eta )D_i{x}_{\max }=(1+\eta )x_i,$$

(29)

where $D_i$ is the dilution factor and $x_i$ is the true number of TSC per reaction for standard i. This shows that every standard will contain approximately the same proportional error as the maximum standard, so the error is systematic and uniform. Now, because Eq. (21) was fitted to the values of ${log}_{10}\left({\tilde{x}}_i\right)$ instead of ${log}_{10}\left(x_i\right)$, the estimated values for field samples produced by Eq. (27) represent ${\tilde{x}}_0=(1+\eta )x_0$ instead of the true $x_0$. Therefore, the estimated number of TSC per reaction in all field samples contains approximately the same proportional error as the estimated number of TSC in the standards. In other words,

$${\displaystyle \frac{{\tilde{x}}_0}{x_0}}={\displaystyle \frac{{\tilde{x}}_i}{x_i}}=1+\eta .$$

(30)

Clearly, it is very important to minimize the error in estimating the number of TSC per reaction in standards, as this error will be passed on unaltered to all estimates for field samples.

To address the other source of error in estimates of the number of TSC per reaction in field samples, we may supplement these estimate with 95% confidence intervals. If assessment of the fitted standard curve and residuals indicates that the assumptions of ordinary least-squares regression are tenable (in particular, normality and homoskedasticity of residuals; see Section 2.4), then a 95% confidence interval for $x_0$ can be estimated using the classical statistical theory of linear calibration (e.g., [30]). The residuals in Figure 4 (left) show no obvious violation of the regression assumptions, so we used these data to calculate examples of 95% confidence intervals for several predicted ${log}_{10}$-transformed initial concentrations calculated from the inverted fitted standard curve. The new $c_{\mathrm{T}}$ values used to predict the ${log}_{10}$-transformed initial concentrations were (arbitrarily) chosen as the means of the observed $c_{\mathrm{T}}$ values for the five different standards plotted in Figure 4 (left). The resulting predicted values and 95% lower and upper confidence limits are displayed in Table 2 and plotted in Figure 4 (right). Note that the actual standard concentrations (vertical dotted lines) fall within the corresponding 95% confidence intervals for the predicted concentrations (gray bars on the horizontal axis). In cases where the regression residuals are approximately normal but clearly heteroskedastic, as is often the case in environmental applications, the classical calibration confidence interval will be only a rough approximation.

The estimated initial copy number given by Eq. (27) represents the number of TSC per reaction in a single well of a multi-well plate. The reaction volume is the volume of S+AM in a single well, but $x_0$ represents simply the number of TSC rather than the number per unit volume. In many applications, it is desirable to convert this value to a concentration representing the number of TSC per unit volume of original sample or standard. To do so, we first divide $x_0$ by reaction volume $V_{\mathrm{r}}$ to convert it to the corresponding concentration in S+AM. We then multiply this concentration by a dimensionless dilution correction factor F that appropriately adjusts it upward for each sample dilution step and downward for each sample concentration step (if any) that occurred during sample processing and analysis, taking care to use the same volume unit (typically either mL or $\mu$L) in all calculations. As a simple example, if the original field sample or standard consitutes 1/5 of the S+AM in each well of a multi-well plate, then the TSC concentration in the original sample or standard would be 5 times the concentration measured in the S+AM, implying $F=5$.

We can estimate amplification factor $\lambda$ using Eq. (23). A simple algebraic manipulation yields

$$\lambda ={10}^{-1/{\beta }_1}={10}^{1/|{\beta }_1|}$$

(31)

[cf. 18, p. 236]. Using this equation and Eq. (3), the proportional amplification efficiency $\epsilon$ can be calculated:

$$\epsilon =\lambda -1={10}^{1/|{\beta }_1|}-1.$$

(32)

This estimate can be used to assess or compare the efficiency of the polymerase chain reaction in different runs (e.g., for quality assurance purposes). If desired, one can also estimate the value of fluorescence ratio $\rho /\phi$ using Eq. (22); the resulting estimates can be used to assess the consistency of the relative values of these two fluorescence parameters in different runs.

A property of Eq. (27) that is important in applications is that $x_0$ depends strongly on fluorescence ratio $\rho /\phi$ and amplification factor $\lambda$. For this reason, qPCR-based estimates of $x_0$ are sensitive to properties of samples that interfere with fluorescence or reduce the efficiency of amplification. Noting that the dependence of $x_0$ on $\rho /\phi$ is linear while its dependence on $\lambda$ is geometric in $c_{\mathrm{T}}$, and recalling that $c_{\mathrm{T}}$ values of roughly 20–35 cycles are typical, it is evident that $x_0$ is particularly sensitive to sample properties that inhibit amplification. This conclusion is strengthened by the fact that sample properties which interfere with fluorescence are likely to affect $\rho$ and $\phi$ similarly, reducing their effect on the ratio.

2.3.4. Accounting for and Minimizing Sample Interference

A variety of potential constituents of environmental samples may interfere with one or more steps of the PCR process, reporter binding, or fluorescence [31,32], typically resulting in an artifactual lowering of the amplification curve, increase in $c_{\mathrm{T}}$, and underestimate of $x_0$. Common methods for reducing interference include sample dilution and use of a specialized master mix designed to reduce inhibition (e.g., TaqMan Environmental Master Mix 2.$0^{\circledR }$, Applied Biosystems, Foster City, California). Additionally, EPA Draft Method C describes a supplemental method for adjusting sample $c_{\mathrm{T}}$ values downward to reduce the effect of inhibition by using purified testes DNA from chum salmon (Oncorhynchus keta) as both a sample processing control and an external positive control [25,26].

2.4. Fitting the Regression Model to Calibration Data

Estimates of intercept and slope parameters ${\beta }_0$ and ${\beta }_1$ are obtained by fitting Eq. (21) to calibration data produced by measuring $c_{\mathrm{T}}$ values for a set of different standards (calibrants), each with several replicates. Standards are solutions containing the target sequence at different concentrations that are prepared or measured in such a way that it is deemed reasonable to treat them as known. A potentially important source of systematic error in qPCR calibration is error in the TSC concentrations assigned to different standards. In practice, these concentrations are often estimated with ddPCR (e.g., EPA Draft Method C [33]). As we explain below in Section 3.2.3, they are therefore dependent on proper determination of the mean droplet volume for the particular laboratory instruments and assay mix employed by the lab performing the ddPCR analysis. In our experience, this determination is rarely done; labs typically accept the default estimate of the mean droplet volume encoded in software supplied by the instrument manufacturer, which is known to be unreliable [34,35,36]. Since ddPCR concentration estimates are inversely proportional to the value of the mean droplet volume (Section 3.2.1), error in the estimated mean droplet volume can result in meaningful error in the resulting estimate of TSC concentration in a standard.

In addition to the problem of accurate estimation of TSC standard concentrations, several purely statistical problems are encountered in fitting a linear regression model to the calibration data. For the estimates of ${\beta }_0$ (intercept) and ${\beta }_1$ (slope) to be valid, it is important to perform various quality control checks on the calibration data to determine if they exhibit excessive variability or a clearly nonlinear trend, or if a small number outliers are present that must be removed. If one wishes not only to obtain point estimates of ${\beta }_0$ and ${\beta }_1$ but also to estimate confidence intervals for them or test statistical hypotheses, then residuals of the fitted model must be examined carefully to assess the key assumptions of the type of regression analysis employed. This procedure often reveals that a different regression model is required, as briefly discussed later in this section.

The calibration data consist of a set of pairs $(u_i,y_i)$ for $i=1,2,3,\dots ,I$, where $u_i$ is the ${log}_{10}$-transformed value of the initial TSC concentration $x_i\left(0\right)$ in the i-th sample (that is, $u_i={log}_{10}\left(x_i\left(0\right)\right)$), $y_i$ is the corresponding observed value of $c_{\mathrm{T}}$ in sample i, and I is the total number of samples (“samples” here are replicates of the different standards, whose concentrations are assumed to be accurately known). In the typical case with equal replication, I is the number of different standards times the number of replicates per standard.

In classical linear regression, one assumes that the observed value of the response variable is a random variable that is the sum of two components: a predicted or “true” value that is a deterministic linear function of the predictor variable (in our case, $u_i$), and a random variable ${\xi }_i$ with mean zero that accounts for the discrepancy between the predicted and observed values. Let $Y_i$ be a random variable representing the observed value of the threshold cycle in the i-th sample, and let ${\widehat{c}}_{\mathrm{T}i}$ denote the corresponding predicted value. Using Eq. (21), ${\widehat{c}}_{\mathrm{T}i}$ is given by

$${\widehat{c}}_{\mathrm{T}i}={\beta }_0+{\beta }_1{u}_i.$$

(33)

The random variable $Y_i$, whose value $y_i$ is observed when $c_{\mathrm{T}}$ is measured in sample i, is the sum of predicted value ${\widehat{c}}_{\mathrm{T}i}$ and random error ${\xi }_i$. That is,

$$\begin{array}{lll}\hfill Y_i& =\hfill & {\widehat{c}}_{\mathrm{T}i}+{\xi }_i\\ & =\hfill & {\beta }_0+{\beta }_1{u}_i+{\xi }_i,i=1,2,3,\dots ,I.\end{array}$$

(34)

We assume the ${\xi }_i$ are independent and have mean zero and finite variance, but we leave open the possibility that the variance of ${\xi }_i$ may be different for different standards. The set of I equations in Eq. (34), together with assumptions regarding the mean and variance of ${\xi }_i$, is the statistical model of the data.

To explain very briefly the process by which Eq. (34) is used to estimate ${\beta }_0$ and ${\beta }_1$, we will consider two different approaches, depending on whether one merely wishes to obtain valid estimates of these parameters and then use them to estimate initial copy number $x_0$ or ${log}_{10}\left(x_0\right)$ (we will call this the simple estimation approach) or also wishes to estimate confidence intervals for the parameters or test statistical hypotheses (we will call this the statistically rigorous approach). Both approaches use the least-squares estimation framework and produce the same estimates, but the simple estimation approach makes no use of probability theory and requires no strong assumptions about statistical properties of the errors. We will begin with the simple estimation approach, then indicate the additional assumptions that must be made and verified if one wishes to employ the statistically rigorous approach.

Both approaches use the quantity ${\beta }_0+{\beta }_1{u}_i$ as the predictor of response variable $Y_i$ for each value $u_i$ of the predictor variable. The difference between the observed value $y_i$ of response variable $Y_i$ and the corresponding predicted value ${\beta }_0+{\beta }_1{u}_i$ is used as a measure of the error contained in the observed value $y_i$ (not the predicted value). To avoid confusing this type of error with random error ${\xi }_i$, it is usually called the residual error, or simply the residual. Thus, the residuals for observations $i=1,2,3,\dots ,I$ can be expressed as

$${\mathrm{residual}}_i={\mathrm{observed}}_i-{\mathrm{predicted}}_i=y_i-({\beta }_0+{\beta }_1{u}_i).$$

(35)

The residual errors have two components: one due to errors in choosing values for ${\beta }_0$ and ${\beta }_1$, and the other due to random errors ${\xi }_i$, which are not observable.

In the ordinary least-squares (OLS) estimation framework, the goal is to find choices of ${\beta }_0$ and ${\beta }_1$ that minimize the sum $Q({\beta }_0,{\beta }_1)$ of squared residuals. For any values of ${\beta }_0$ and ${\beta }_1$, this sum is given by

$$Q({\beta }_0,{\beta }_1)=\sum _{i=1}^I{\left({\mathrm{residual}}_i\right)}^2=\sum _{i=1}^I{[y_i-({\beta }_0+{\beta }_1{u}_i)]}^2.$$

(36)

Least-squares estimates of the intercept and slope parameters, denoted by ${\widehat{\beta }}_0$ and ${\widehat{\beta }}_1$, are simply the values that minimize $Q({\beta }_0,{\beta }_1)$ and are easily found. The objective, then, is to minimize the contribution that errors in the values of ${\beta }_0$ and ${\beta }_1$ make to the total residual error, leaving random errors ${\xi }_i$ as the dominant contributors. The residuals should then be reasonable estimates of random errors ${\xi }_i$ if the underlying relationship between predictor and response variables truly is linear. With the simple estimation approach, it is unnecessary to make assumptions about statistical properties of the errors. In particular, we need not make assumptions about the specific probability distribution they come from, whether their variance is homogeneous for all $u_i$, or even whether they are independent. This approach, then, is simply a “curve-fitting” approach.

The statistically rigorous approach to estimating ${\beta }_0$ and ${\beta }_1$ finds the estimates in exactly the same way but focuses on random variables $Y_i$ and ${\xi }_i$ and the probability distributions they are assumed to come from. It imposes several specializing assumptions about statistical properties of errors ${\xi }_i$ that make it possible to estimate confidence intervals for ${\beta }_0$ and ${\beta }_1$ and, if desired, test statistical hypotheses about them. Specifically, the statistically rigorous approach assumes the errors ${\xi }_i$ are independent and (approximately) normally distributed with constant mean zero and with a variance that is positive and finite for all i, constant for all i corresponding to any given standard, but possibly different for different standards. If the qPCR analysis was performed correctly, the trend in the data should be well described by a straight line (as in Figure 4 above), in which case the residuals should show no compelling evidence of serial correlation (which would violate the independence assumption). The normality assumption should be checked visually with a normal quantile-quantile plot (optionally supplemented with a formal test, such as the Shapiro-Wilk test) and in our experience usually is tenable. The simplest assumption regarding the error variance is that it is homogeneous for all standards, in which case the variance of the residuals should be homogeneous, as well. This property can be checked visually with a plot of residuals versus ${log}_{10}$-transformed standard concentration and a box plot (optionally supplemented with a Levene test for homogeneity) and often is dubious or clearly untenable.

Unless there is strong evidence to the contrary, the variance homogeneity assumption should be accepted as valid and OLS regression should be employed. If, however, variance is clearly not homogeneous across standards, weighted least-squares (WLS) regression is sometimes recommended (this is, in fact, the default method used in the spreadsheet workbook that accompanies EPA Draft Method C). With this type of regression, the squared residuals are multiplied by positive weights before summing. The weighted sum $Q({\beta }_0,{\beta }_1)$ of squared residuals is then given by

$$Q({\beta }_0,{\beta }_1)=\sum _{i=1}^I{w}_i·{\left({\mathrm{residual}}_i\right)}^2=\sum _{i=1}^I{w}_i{[y_i-({\beta }_0+{\beta }_1{u}_i)]}^2$$

(37)

(e.g., [37], Section 2.1.3.1), where $w_i>0$ is the weight attached to the i-th squared residual. With the weights fixed, least-squares estimates of ${\beta }_0$ and ${\beta }_1$ are found by minimizing $Q({\beta }_0,{\beta }_1)$, as with OLS regression.

Moving $w_i$ inside the brackets in Eq. (37), $Q({\beta }_0,{\beta }_1)$ can be written in the equivalent form,

$$Q({\beta }_0,{\beta }_1)=\sum _{i=1}^I{[{\tilde{y}}_i-({\beta }_0\sqrt{w_i}+{\beta }_1{\tilde{u}}_i)]}^2,$$

(38)

where

$${\tilde{y}}_i=y_i\sqrt{w_i},{\tilde{u}}_i=u_i\sqrt{w_i}.$$

(39)

If we view ${\tilde{y}}_i$ as a response variable and $\sqrt{w_i}$ and ${\tilde{u}}_i$ as predictor variables, then Eq. (38) is just the sum of squared residuals for a two-predictor OLS regression model with response variable ${\tilde{y}}_i$, predictor variables $\sqrt{w_i}$ and ${\tilde{u}}_i$, slope parameters ${\beta }_0$ and ${\beta }_1$, and fixed intercept 0. The underlying model of the data is

$${\tilde{Y}}_i={\beta }_0\sqrt{w_i}+{\beta }_1{\tilde{u}}_i+\sqrt{w_i}{\xi }_i,$$

(40)

where the last term on the right is the error term. This 2-predictor model is used to estimate ${\beta }_0$ and ${\beta }_1$ and their confidence intervals. The variance of the error term is $w_i{\sigma }_i^2$, which can be made approximately constant for all i if the weights are chosen so that $w_i$ is inversely proportional to ${\sigma }_i^2$.

The residuals of the 2-predictor model for observations $1,2,3,\dots ,n$ are given by

$${\mathrm{residual}}_i={\tilde{y}}_i-({\beta }_0\sqrt{w_i}+{\beta }_1{\tilde{u}}_i).$$

(41)

Once parameters ${\beta }_0$ and ${\beta }_1$ have been estimated, these residuals should be reasonable estimates of errors $\sqrt{w_i}{\xi }_i$ if the weights are chosen well, and they therefore should be approximately homogeneous over i. If so, valid confidence intervals for ${\beta }_0$ and ${\beta }_1$ can be estimated. However, because the variance of ${\xi }_i$ is often relatively low for high standards and high for low standards, greater weight will be attached to squared residuals for high standards that for low ones, potentially causing the fitted least-squares regression model to fit high standards better than low ones. As result, the accuracy of predicted target sequence copy numbers in samples with low concentrations may be lower than the accuracy at high concentrations (an example with real data is given below). In that case, if one is mainly interested in accurately estimating target sequence copy numbers in environmental samples and if some of these are likely to be low, then use of OLS regression may be preferable.

2.5. Example: Comparing Statistical Properties of fitted OLS and WLS Standard Curves

A fundamental principle of qPCR analysis is that measuring target-sequence copy numbers in environmental samples depends critically on a standard curve that has been fitted to a set of calibration data comprising pairs $(u_i,y_i)$ of known ${log}_{10}$-transformed initial copy numbers $u_i={log}_{10}\left(x_i\left(0\right)\right)$ and corresponding measured threshold cycles $y_i$ for a set of standards with different concentrations of the target sequence. The fitted standard curve is linear and is therefore determined by its intercept and slope parameters ${\beta }_0$ and ${\beta }_1$. The least-squares estimates of these parameters, when inserted in the inverted standard curve, determine the estimated target-sequence copy number for each $c_{\mathrm{T}}$ value measured in an environmental sample. The purely statistical issue of how best to fit a standard curve to calibration data to obtain reasonable estimates of ${\beta }_0$ and ${\beta }_1$ is therefore of great importance in qPCR analysis but is not as simple as one might hope.

EPA Draft Method C has been used extensively in the United States for monitoring $E.coli$ contamination of aquatic habitats. A spreadsheet workbook distributed with this method performs various calculations that are required for assessing data quality, fitting standard curves, and calculating copy numbers in environmental samples. To simplify use of the workbook by laboratory personnel with little statistical background, the current version of this workbook automatically uses WLS regression to fit standard curves to calibration data, without basing this choice of regression approach on a careful assessment of the assumptions of OLS and WLS regression or the resulting model fits. Additionally, the spreadsheet automatically employs as weights the known ${log}_{10}$-transformed copy numbers for different standards, with no assessment of whether these weights successfully homogenize the variance of residuals for the fitted WLS model. The rationale for this choice of weights is that one often finds the variance of residuals to be higher for the lower standard concentrations than for the higher concentrations, so it is reasonable to think that standard concentrations might be approximately inversely proportional to the variance of the corresponding residuals. But like any other choice of weights, this choice requires careful assessment to determine whether it succeeds in homogenizing the variance of residuals and whether the resulting intercept and slope parameters perform better than (or as well as) those for standard OLS regression when used in Eq. (26) or (27) for estimating copy numbers in samples.

We now consider an example with real qPCR calibration data that allows us to compare the performance of OLS and two types of WLS regression which use different weighting schemes, including the scheme used in the current EPA Method C workbook. The results of our comparison are specific to the particular data set we employ but illustrate the importance of assessing alternative regression models and their assumptions when fitting qPCR standard curves.

The data set we consider consists of 6 runs of a qPCR calibration study. Each run used 5 standards (i.e., 5 different known concentrations of a particular DNA target sequence), with 3 replicates of each. In each run, the threshold cycle $c_{\mathrm{T}}$ was determined for each replicate of each standard, yielding 6 separate sets of calibration data. The spreadsheet workbook automatically runs preliminary quality-assurance checks to determine (among other things) whether the data sets for different runs are sufficiently consistent to be combined. In this example, the spreadsheet concluded that it was appropriate to combine the data from all 6 runs, yielding a composite data set with 18 replicates (6 runs × 3 replicates per run) for each of the 5 standards. The workbook then used this composite data set to estimate intercept and slope parameters ${\beta }_0$ and ${\beta }_1$ and their 95% confidence intervals. EPA Draft Method C refers to the resulting standard curve as a “composite curve”.

We assessed three procedures for estimating ${\beta }_0$ and ${\beta }_1$ with these data:

• Simple OLS regression
• WLS regression using the EPA Draft Method C weights
• WLS regression using alternative weights based on the data.

We will refer to these procedures as OLS, EPA WLS, and Alternative WLS regression. As mentioned above, the EPA WLS weights are the known ${log}_{10}$-transformed copy numbers for the different standards. The Alternative WLS weights are $1/s_j^2$, where $s_j$ is the sample standard deviation of the 18 residuals for standard j in the OLS model or, equivalently, of the 18 measured $c_{\mathrm{T}}$ values (e.g., [37], Section 2.1.3.1). We note that our parameter estimates and 95% confidence intervals for EPA WLS regression, computed with our own program written in the R programming language [38], agree exactly with results calculated by the EPA Draft Method C spreadsheet workbook.

The assumptions of normality and variance homogeneity of residuals were assessed both visually and with formal statistical tests. Normality was assessed using quantile-quantile plots, supplemented with Shapiro-Wilk tests. Variance homogeneity was assessed using box plots, supplemented with Levene tests. The overall fit of the different regression models was quantified on the original measurement scale using the mean absolute error MAE, defined by

$$\mathrm{MAE}={\displaystyle \frac{1}{n}}\sum _{i=1}^I\left|{\mathrm{residual}}_i\right|,$$

which avoids squaring the residuals and is therefore less sensitive to outliers than is the root mean squared error.

Figure 5 shows the results of our assessment. The first column of panels applies to OLS regression, the second the EPA WLS regression, and the third to Alternative WLS regression. The first row of panels shows $c_{\mathrm{T}}$ versus ${log}_{10}\left(x_0\right)$ calibration data and fitted regression lines, the second shows residuals as defined by Eq. (35) for OLS regression and by Eq. (41) for WLS regression, the third shows corresponding box plots of the residuals, the fourth shows normal quantile-quantile plots of the residuals, and the last shows inverse predictions of ${log}_{10}\left(x_0\right)$ for the five standards, based on means of the corresponding measured $c_{\mathrm{T}}$ values and Eq. (24).

Several interesting conclusions emerge from these results. The first row of panels suggest that all three regression models fit the calibration data well at high standard concentrations, but that the Alternative WLS model does not fit the data at the two lowest concentrations as well as the OLS and EPA WLS models do. The second and third rows of panels show that the variance homogeneity assumption for residuals is not tenable for either the OLS or EPA WLS models but clearly is for the Alternative WLS model. The first panel in the second row also shows that the standard deviations $s_j$ of $c_{\mathrm{T}}$ values listed at the bottom of the panel (these are also the standard deviations of residuals for the OLS model) are smallest for the middle and highest standards and much larger for the two lowest standards. Because the Alternative WLS model weights residuals for each standard by $1/s_j^2$, residuals for the middle and highest standards are weighted most heavily and those for the two lowest standards are weighted least heavily. This is why the Alternative WLS model fits the middle and highest standards particularly well and the lowest standard particularly poorly. The fourth row of panels shows that the normality assumption is tenable for residuals of all three regression models. In the last row of panels, the means of the measured $c_{\mathrm{T}}$ values for the different standards were inserted in the inverse prediction equation stated in Eq. (24) to predict ${log}_{10}$ copy numbers for the standards. All of these predicted ${log}_{10}$ copy numbers should be close to the corresponding actual ${log}_{10}$ copy numbers in the standards if the calibration data are of good quality and the standard curve has been fitted well. Note that the OLS and EPA WLS models predict well at the lowest and two highest concentrations but noticeably worse at the two intermediate concentrations. The Alternative WLS model predicts well at the three highest concentrations but noticeably worse at the two lowest concentrations. These patterns in row 5 are consequences of the patterns in row 1 and result from the choice of regression model and weights used to fit the standard curve.

We conclude that, in this particular example using real data, it is not clear that using either type of WLS regression model is preferable to using the OLS model. The OLS and EPA WLS models performed similarly with regard to predicting TSC concentrations across the entire range of standards, but in neither case was the variance of residuals homogeneous, as required for proper parametric estimation of confidence intervals or tests of statistical hypotheses. By contrast, the Alternative WLS regression model successfully homogenized the variance of residuals but predicted concentrations less accurately than the OLS or EPA WLS models at low concentrations. When such cases arise in environmental applications, it seems more important to obtain accurate estimates of sample TSC concentrations across the entire range of standards at the cost of some heterogeneity of the regression residuals, than it is to resolve the problem of variance heterogeneity at the cost of inaccurate estimates for low TSC concentrations.

These conclusions, of course, apply only to the particular example considered. More generally, we may tentatively conclude that neither the desirability of using WLS regression nor the best weighting scheme to employ should simply be assumed; both choices should be assessed on a case-by-case basis. Because of the importance of proper fitting of the standard curve in qPCR analysis, we think a more thorough statistical assessment of alternative fitting methods is needed. Such an assessment, however, is beyond the scope of the present study.

3. Droplet Digital PCR

3.1. A Brief Overview of ddPCR

As in the case of qPCR, numerous protocols are available for ddPCR analysis of environmental samples. (An example of a ddPCR workflow that has been widely used in the state of Michigan for monitoring SARS-CoV-2 concentrations in wastewater is presented in Appendix C.) Regardless of the protocol used, the first step in ddPCR analysis is sample preparation, which releases the genetic material contained in cells of the target taxon and cleans the sample to varying degrees for subsequent analysis. The methods used are the same as those described for qPCR in Section 2.1 and will not be repeated here.

Figure 5. Calibration results for a composite data set that includes data from 6 qPCR calibration runs, each with 3 replicates of the same 5 standards (total replicates per standard = 18). Standard concentrations are the same as in Figure 4. Three regression models were assessed and are listed at the top of the figure. Row 1: Data and fitted standard curves, with estimates of the intercept, slope, and mean absolute error (MAE). Row 2: Plots of regression residuals versus ${log}_{10}\left(x_0\right)$. Residuals for WLS models are those for the equivalent 2-predictor model in Eq. (40), since this is the model used to estimate parameter values and confidence intervals. Standard deviations of the residuals for each standard are listed at the bottom of the panel. Row 3: Box plots of regression residuals versus ${log}_{10}\left(x_0\right)$. P-values for Levene’s test are also shown. Row 4: Normal quantile-quantile plots of residuals, with 95% confidence envelopes. P-values for the Shapiro-Wilk test are also shown. Row 5: Inverse prediction of ${log}_{10}$ initial copy numbers of the five standards (filled yellow circles on the horizontal axis), based on means of the corresponding measured $c_{\mathrm{T}}$ values (filled yellow circes on the vertical axis). Dotted vertical lines: actual ${log}_{10}$ initial copy numbers of the standards. Data: John J. Hart, Robert B. Annis Water Resources Institute.

Figure 5. Calibration results for a composite data set that includes data from 6 qPCR calibration runs, each with 3 replicates of the same 5 standards (total replicates per standard = 18). Standard concentrations are the same as in Figure 4. Three regression models were assessed and are listed at the top of the figure. Row 1: Data and fitted standard curves, with estimates of the intercept, slope, and mean absolute error (MAE). Row 2: Plots of regression residuals versus ${log}_{10}\left(x_0\right)$. Residuals for WLS models are those for the equivalent 2-predictor model in Eq. (40), since this is the model used to estimate parameter values and confidence intervals. Standard deviations of the residuals for each standard are listed at the bottom of the panel. Row 3: Box plots of regression residuals versus ${log}_{10}\left(x_0\right)$. P-values for Levene’s test are also shown. Row 4: Normal quantile-quantile plots of residuals, with 95% confidence envelopes. P-values for the Shapiro-Wilk test are also shown. Row 5: Inverse prediction of ${log}_{10}$ initial copy numbers of the five standards (filled yellow circles on the horizontal axis), based on means of the corresponding measured $c_{\mathrm{T}}$ values (filled yellow circes on the vertical axis). Dotted vertical lines: actual ${log}_{10}$ initial copy numbers of the standards. Data: John J. Hart, Robert B. Annis Water Resources Institute.

Droplet digital PCR (ddPCR) is an endpoint form of PCR analysis where each sample is partitioned into roughly 10,000–20,000 nanoliter-size droplets, each functioning as a separate PCR reactor. As in qPCR, the prepared sample is combined with an assay mix containing primers, probes, dNTPs, DNA polymerase, and other constituents, with primers and probes binding specifically to the target sequence of DNA. But unlike qPCR, fluorescence of droplets is not monitored in real time; it is assessed only once, after the number of amplification cycles is large enough so droplets that fluoresce above background can be reliably distinguished from those that do not (typically 40–60 cycles). An individual droplet will fluoresce above background if and only if it originally contained one or more copies of the target DNA sequence (which was then amplified). As will be shown in Section 3.2, the original concentration of the target sequence can be estimated from two simple pieces of information: the observed proportion of droplets that did not fluoresce above background and the mean droplet volume.

In order to understand how ddPCR works, it is essential to have a basic idea of the manner in which droplets of the mixture of prepared sample and ddPCR assay mix are created. Some of the details vary between different ddPCR systems, but the process is as follows for the ${\mathrm{Bio-Rad}}^{\circledR }{\mathrm{QX200}}^{\mathrm{TM}}$ PCR system. Very briefly, an instrument called a droplet generator uses pressurized air to push 20 $\mu$L of aqueous S+AM and 60 $\mu$L of a special oil through three separate microchannels (width and depth typically less than 100 $\mu$m) in a plastic chip. Oil flows through two of the microchannels and S+AM through the third. The two oil microchannels join the S+AM microchannel from opposite sides at a common point (Figure 6). As the flowing oil and S+AM merge, the fact that the two liquids are immiscible becomes important: instead of mixing, droplets of S+AM are sequentially pinched off by oil flowing in from opposite sides of the S+AM microchannel, then remain suspended within the oil as it flows through an outflow microchannel. This process partitions the 20 $\mu$L of S+AM into roughly 20,000 nanoliter-size droplets, each with a volume of approximately 1.0 nL (the actual volume is somewhat variable but typically is closer to 0.8 nL [34,35,36]).

When nearly the entire 20 $\mu$L of S+AM has been processed by the droplet generator, the emulsion of aqueous droplets in oil is transferred to a multi-well plate and placed on a thermocycler for PCR amplification, which occurs separately in each droplet that contains one or more copies of the target sequence. Once amplification is complete, the plate is placed on an instrument called a droplet reader. The emulsion in each well is aspirated sequentially by a needle and passed through a microchannel where the droplets pass, one at a time, through a beam of light at the excitation wavelength for the reporter being employed. Fluourescence of the droplet-in-oil emulsion at the emission wavelength of the reporter is continuously monitored. Software for the droplet reader employs an automated peak-detection algorithm that assesses each putative fluorescence peak to determine whether to accept it as corresponding to a valid droplet, and the fluorescence amplitude of each accepted peak is recorded.

Even droplets that contain no TSC may have a positive but small fluorescence amplitude, so the analyst must set a threshold level that reliably separates flourescence amplitudes in positive control droplets from fluorescence amplitudes in negative control and no-template control droplets. Software accompanying the droplet reader uses this threshold to classify each S+AM droplet as either positive (fluorescing above the threshold) or negative (not fluorescing above the threshold) for the target sequence. The software then uses the observed proportion of negative droplets and an assumed mean droplet volume to calculate an estimate of the target sequence concentration as TSC per 20 $\mu$L of S+AM (as discussed below, the reported concentration will require adjustment if there is a meaningful difference between the mean droplet volume coded into the instrument software and the actual mean volume of droplets produced by the particular droplet generator and master mix used in a given study).

Depending on the purpose of the study, it may be necessary to convert concentrations in the S+AM to concentrations in the original field samples, accounting for all concentration and dilution steps during sample preparation and analysis. Unlike qPCR, ddPCR does not require a standard curve and does not depend on estimates of amplification efficiency. These are important advantages.

3.2. Basic Equations of the ddPCR Process

As discussed in Section 3.1, the concentration of the target sequence in a sample can be estimated with ddPCR using only the proportion of negative droplets, as determined by the droplet reader, and the mean droplet volume. We now show why this is true. Our exposition is somewhat more detailed than the various expositions we have encountered in the literature and on instrument manufacturer and vendor websites, since we feel it is important to be aware of the key assumptions and approximations underlying the theory and the manner in which the main results can be logically derived.

The key result to be justified is the formula typically used in ddPCR analysis to estimate the concentration of TSC in the S+AM (number per volume of S+AM):

$$C=-{log}_e(m_0/m)/\overline{v}$$

(42)

[e.g., 20, pp. 54–55]). Here, C is the TSC concentration in the S+AM, m is the total number of accepted droplets scanned by the droplet reader, $m_0$ is the number of these droplets that were classified as negative (i.e., that did not fluoresce above background and therefore contained no TSC), $m_0/m$ is the proportion of droplets that were negative, and $\overline{v}$ is the mean droplet volume. It is important to be aware that the mean droplet volume depends on the type of droplet generator and master mix used (and possibly on other factors that have not yet been assessed), with most of the reported optically-measured values ranging from 0.7 and 0.9 nL [34,35,36,39].

In our derivation of Eq. (42), we will assume that S+AM, oil, and droplets suspended in oil all move through the various microchannels shown in Figure 6 via plug flow. Thus, each type of fluid is assumed to move as if on a conveyor belt, and any arbitrary segment of the fluid (and any TSC it contains) remains intact as it moves through the microchannel, like a plug or piston. Consistent with the plug flow assumption, it is useful to think of the stream of S+AM flowing through the sample microchannel as consisting of a sequence of disjoint segments, each having approximately the same volume and each becoming a single droplet when pinched off by inflow of oil from the two oil microchannels, as illustrated in Figure 6.

3.2.1. Estimating the TSC Concentration

We begin with several assumptions regarding probabilistic properties of the number of TSC in different segments of the S+AM as it flows through the sample microchannel, prior to the inflow of oil (Figure 6):

• Before being drawn into the sample microchannel, the S+AM is homogeneous (well mixed) in the sense that any given TSC that was present in the original sample is equally likely to be contained in any subset of the S+AM of fixed volume $v\ll V$, where V is the total S+AM volume ($\approx 20\mu$L).
• The droplet generator partitions the entire S+AM volume V (or nearly so) into a large number n of droplets ($n\approx 20,000$ is typically assumed) of variable but similar volume.
• The process by which the N TSC present in the S+AM are allocated to the n disjoint segments of S+AM flowing through the sample microchannel, and hence to the n resulting droplets, is a stochastic partition process equivalent to independent random allocation of N objects (the TSC) to n boxes (the segments or droplets) labeled $1,2,\dots ,n$, where each object has probability ${\pi }_j$ of being allocated to box j (Figure 7).
• The probability ${\pi }_j$ that any given TSC is allocated to S+AM segment or droplet j is simply the fraction of S+AM volume V that the segment or droplet contains, so that

  $${\pi }_j=v_j/V,j=1,2,\dots ,n.$$

  (43)
• In order to arrive at the simple standard formula for TSC concentration C in the S+AM stated in Eq. (42), we must make the additional assumption that the volumes of different droplets are similar enough so the approximation,

  $$v_j\approx \overline{v},j=1,2,\dots ,n,$$

  (44)

  is adequate, in which case it follows that

  $${\pi }_j\approx \overline{v}/V:=\overline{\pi },j=1,2,\dots ,n.$$

  (45)

These assumptions imply that the numbers $Z_1,Z_2,\dots ,Z_n$ of TSC in droplets (and segments) $1,2,\dots ,n$ are random variables that jointly have a multinomial distribution $M_n(N;{\pi }_1,{\pi }_2,\dots ,{\pi }_n)$ with probability function

$$\begin{array}{}\mathrm{(46)}& \hfill Pr(Z_1=z_1,Z_2=z_2,\dots ,Z_n=z_n)& =\hfill & {\displaystyle \frac{N!}{z_1!z_2!\cdots z_n!}}{\pi }_1^{z_1}{\pi }_2^{z_2}\cdots {\pi }_n^{z_n}\mathrm{(47)}& & \approx \hfill & {\displaystyle \frac{N!}{z_1!z_2!\cdots z_n!}}{\overline{\pi }}^N,\end{array}$$

where realized values $z_1,z_2,\dots ,z_n$ are constrained to sum to N and we used (in the second step) the assumption that ${\pi }_j\approx \overline{\pi }$ for every droplet j. Symbols and associated dimensions employed in the qPCR equations here and below are summarized in Table 3.

An important property of the multinomial distribution (and the partition process it represents here) is that the numbers of objects allocated to different boxes are not independent, because the numbers of objects and boxes are fixed and finite (e.g., in the extreme case where all N objects happen to be allocated to a particular box $j^{\prime }$, none can be allocated to any other box). The correlation between the numbers of TSC in any two S+AM droplets is negative and is given by

$$\mathrm{Cor}(X_i,X_j)=-\sqrt{{\displaystyle \frac{\phantom{t}{\pi }_i{\pi }_j\phantom{t}}{(1-{\pi }_i)(1-{\pi }_j)}}},i\ne j$$

(48)

(e.g., [40]). Under our assumptions, then,

$$\mathrm{Cor}(X_i,X_j)\approx -\overline{\pi }/(1-\overline{\pi })=-1/(n-1).$$

(49)

Thus, for $n={10}^k$ with $k\ge 3$, we find that $\mathrm{Cor}(X_i,X_j)\approx -{10}^{-k}$. It follows that the between-droplet correlation in TSC counts is negligible for $n>1000$, which clearly holds for typical ddPCR applications with $n\approx 20,000$.

The key implication of this result is that, because the number n of droplets is so large, the number of TSC in different droplets can be assumed to be independent random variables, even if the entire 20 $\mu$L of S+AM is partitioned to produce the droplets. Under the independence assumption, the number of TSC in each droplet when it was created by the droplet generator has the same probability distribution and is not influenced by the number of TSC in any other droplet. This distribution is easily shown to be binomial with “success” probability $\pi =1/n$ and number of “trials” N (hence, the mean number of TSC per droplet = $N\pi =N/n$), but as we will see shortly, it is only necessary to determine the probability that any given S+AM droplet contains no TSC; the rest of the binomial distribution plays no role in estimating the concentration of TSC, nor do we have any way of knowing the number of TSC that were present in individual positive droplets when they were created.

The probability $P_0$ that any given S+AM droplet contains no TSC is simply the probability that none of the N TSC was allocated to the segment from which the droplet was formed. Under the independence assumption, and using Eq. (43), $P_0$ is given by

$$P_0={(1-\pi )}^N={(1-1/n)}^N.$$

(50)

Taking the logarithm of both sides and solving for N yields the following formula for the total number of TSC in the S+AM:

$$N={\displaystyle \frac{{log}_e\left(P_0\right)}{{log}_e(1-1/n)}}.$$

(51)

This formula actually holds for logarithms with respect to any base (provided the base is the same in the numerator and denominator), but derivation of Eq. (42) is simpler with base e.

The concentration C of TSC in the S+AM is the number N of TSC divided by the total volume V of S+AM. That is,

$$C={\displaystyle \frac{N}{V}}={\displaystyle \frac{{log}_e\left(P_0\right)}{V{log}_e(1-1/n)}}={\displaystyle \frac{{log}_e\left(P_0\right)}{\left({\sum }_{i=1}^n{v}_i\right){log}_e(1-1/n)}}.$$

(52)

Because n is so large (≈ 20,000), $1/n$ is very small (≈ 0.00005). To a very good approximation, then, we may assume

$${log}_e(1-1/n)=-1/n$$

(53)

(MacLaurin series with remainder). Eq. (52) then becomes

$$C={\displaystyle \frac{-{log}_e\left(P_0\right)}{(1/n){\sum }_{i=1}^n{v}_i}}={\displaystyle \frac{-{log}_e\left(P_0\right)}{\overline{v}}}.$$

(54)

In practice, not all of the n droplets into which the S+AM was partitioned are scanned and accepted by the droplet reader. Therefore, let $m<n$ denote the number of droplets that are scanned and accepted. To obtain a valid estimate of the TSC concentration in the S+AM with this subset of droplets, the m droplets in the subset must be a representative sample of the total number n of droplets in the S+AM. Because both m and n typically are large ($>{10}^4$) and there is no apparent source of bias in determining which droplets will be read, it is reasonable to assume this requirement is satisfied. The following corollary properties then hold:

• The probability $P_0^{\prime }$ that a randomly chosen droplet from the m scanned and accepted droplets contains no TSC is the same as the probability $P_0$ that a randomly chosen droplet from the full set of n droplets contains no TSC. That is, $P_0^{\prime }=P_0.$
• The ratio $N^{\prime }/m$ of TSC to droplets in the m scanned and accepted droplets is the same as the ratio $N/n$ in the full set of n droplets. That is, $N^{\prime }/m=N/n$.
• The mean droplet volume ${\overline{v}}^{\prime }$ in the m scanned and accepted droplets is the same as the mean droplet volume $\overline{v}$ in the full set of n droplets. That is, ${\overline{v}}^{\prime }=\overline{v}$.
• It follows from properties 1 and 2 that, to a very good approximation, the ratio $N^{\prime }/V^{\prime }$ of the cumulative number $N^{\prime }$ of TSC in the m scanned and accepted droplets to the cumulative volume $V^{\prime }$ of those droplets is the same as the ratio $N/V$ of the cumulative number N of TSC in the total number n of droplets (and in the S+AM)to the total volume V of those droplets (and of the S+AM). That is,

  $$C={\displaystyle \frac{N}{V}}={\displaystyle \frac{N^{\prime }}{V^{\prime }}}=C^{\prime },$$

  (55)

  where C is the TSC concentration in the combined n droplets (and in the S+AM) and $C^{\prime }$ is the combined TSC concentration in the m droplets scanned and accepted by the droplet reader.

Properties 1, 3, and 4 imply that we may substitute $P_0^{\prime }$ for $P_0$ and ${\overline{v}}^{\prime }$ for $\overline{v}$ in Eq. (54), and that to a very good approximation, $C=C^{\prime }$. Eq. (54) can then be written as

$$C=C^{\prime }=-{log}_e\left(P_0^{\prime }\right)/{\overline{v}}^{\prime }=-{log}_e\left(P_0^{\prime }\right)/\overline{v}.$$

(56)

In Section 3.2.2 below, we show that the maximum-likelihood estimator ${\widehat{P}}_0^{\prime }$ for probability $P_0^{\prime }$ is $m_0/m$, where $m_0$ is the number of the m droplets scanned and accepted by the droplet reader that did not fluoresce above the threshold. Substituting ${\widehat{P}}_0^{\prime }$ for $P_0^{\prime }$ in Eq. (56), the formula for estimating the concentration C of TSC in the S+AM becomes

$$C=C^{\prime }=-{log}_e(m_0/m)/\overline{v},$$

(57)

which agrees with Eq. (42).

As in the case of qPCR, it is often desirable to convert this concentration in S+AM to the corresponding concentration in the original sample or standard. The process for doing so is the same: we multiply the concentration in S+AM by a dimensionless dilution correction factor F that adjusts it appropriately to account for sample dilution and concentration steps that occurred during sample processing and analysis (for an example, see [34]). The concentration $C_{\mathrm{sample}}$ in the original sample is therefore given by

$$C_{\mathrm{sample}}=FC=-{log}_e\left(P_0^{\prime }\right)F/\overline{v}.$$

(58)

In order to use Eq. (56) to estimate the TSC concentration in the S+AM, the values of parameters $\overline{v}$ and $P_0$ must be determined. The value of $P_0$ can be determined statistically using a maximum-likelihood estimator (Section 3.2.2). An estimate of $\overline{v}$ is coded into software supplied by instrument manufacturers, but such universal estimates have been found to be unreliable. More specifically, it has been shown by multiple laboratories and standards or metrology organizations around the world that the value of $\overline{v}$ differs meaningfully for different droplet generators and master mixes [34,35,36]. Thus, while it is true that ddPCR does not require a standard curve, it does require calibration. The recommended procedure is to estimate $\overline{v}$ for each combination of droplet generator and master mix used in a particular laboratory by measuring droplets under a microscope in accordance with published protocols (e.g., [34,35]). A simpler alternative, which we propose in Section 3.2.3, is to employ one-point calibration using a certified reference material (CRM) with known TSC concentration.

Before discussing estimation of $\overline{v}$ and $P_0$, we perhaps should explain why we have not mentioned what is referred to as “Poisson statistics” in ddPCR documentation on manufacturer websites. In the context of problems like the one we are considering, the standard theory of spatial Poisson stochastic processes does not apply, because TSC are allocated to droplets by partitioning a fixed number of TSC. Instead, the Poisson distribution arises as an approximation to the binomial distribution, which in turn is an approximation to the underlying multinomial distribution. The binomial approximation should be very good, because the extremely large number n of droplets ensures that the correlation between TSC numbers in different droplets is negligible, and variation in droplet size (hence, ${\pi }_j$) is minimal (see Section 3.2.2). Because n is so large, the Poisson approximation also should be fairly good, provided the average number of TSC per droplet ($N/n$) is less than about 10 [41]. However, it is actually only necessary for the Poisson approximation to $P_0^{\prime }$ to closely approximate the corresponding binomial probability given by Eq. (50) with $P_0=P_0^{\prime }$. The Poisson approximation to $P_0^{\prime }$ is

$$P_0^{\prime }=e^{-N^{\prime }/m}=e^{-N/n}.$$

(59)

Numerical analysis of the difference between the binomial probabiity in Eq. (50) and the Poisson approximation to it in Eq. (59) show that for $n\ge 10,000$, the difference is less than 0.00002 for all values of TSC number $N\ge 1$. However, we see no reason to use the Poisson approximation instead of the binomial probability, which is simple to compute and has a more transparent relationship with the underlying multinomial distribution.

3.2.2. Estimating the Probability that a Droplet Has no TSC

We now outline the statistical theory for estimating $P_0^{\prime }$, the probability that any given scanned droplet contains no TSC, based on droplet-reader results for a large number $m<n$ of scanned and accepted droplets. Let $Y_j$ be a binary (Bernoulli) random variable indicating whether droplet j is negative ($Y_j=0$) or positive ($Y_j=1$) when scanned. Thus, for each droplet j, we will have either $Y_j=0$, if the droplet contained no TSC when it was formed, or $Y_j=1$, if the droplet contained one or more TSC when it was formed (which were then amplified by PCR before being scanned by the droplet reader). The probability distribution of $Y_j$ can therefore be expressed as

$$Y_j=\left\{\begin{array}{c}0,\mathrm{with}\mathrm{probability}{P}_0^{\prime }\hfill \\ 1,\mathrm{with}\mathrm{probability}1-P_0^{\prime }\hfill \end{array}\right.$$

(60)

for all j.

Recalling the independence assumption justified by the multinomial argument presented in Section 3.2.1, we may treat $Y_1,Y_2,Y_3,\cdots$ as independent random variables. Without loss of generality, we may interpret $Y_j=0$ as a “success” and $Y_j=1$ as a “failure”. Then by definition of the binomial distribution, the number of negative droplets in m accepted droplets has a binomial distribution with success probability $P_0^{\prime }$ and number of trials m. (Note: the number of positive droplets is also binomial but with success probability $1-P_0^{\prime }$; we see no valid reason to approximate the binomial with a Poisson distribution, as is sometimes done, since this inevitably introduces additional error and all the necessary statistical theory is available for the binomial distribution.) It follows that the likelihood $L\left(P_0^{\prime }\right)$ of obtaining the observed number $m_0$ of successes in m Bernoulli trials ($0\le m_0\le m$) is given by

$$L\left(P_0^{\prime }\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{m}{m_0}}\right)P_0^{\prime m_0}{(1-P_0^{\prime })}^{m-m_0}.$$

(61)

We wish to choose $P_0^{\prime }$ to maximize this likelihood. Note that $L\left(0\right)=0=L\left(1\right)$ and $L\left(P_0^{\prime }\right)>0$ for all $P_0^{\prime }\in (0,1)$, suggesting there is an interior maximum. The log-likelihood function $\ell \left(P_0^{\prime }\right):={log}_e\left(L\left(P_0^{\prime }\right)\right)$ is given by

$$\ell \left(P_0^{\prime }\right)=\mathrm{const}+m_0{log}_e\left(P_0^{\prime }\right)+(m-m_0){log}_e(1-P_0^{\prime }).$$

(62)

Differentiating $\ell \left(P_0^{\prime }\right)$ with respect to $P_0$ and setting the derivative to zero (and confirming the second derivative at the implied value of $P_0^{\prime }$ is negative), we find that the maximum-likelihood estimator ${\widehat{P}}_0^{\prime }$ for $P_0^{\prime }$ is

$${\widehat{P}}_0^{\prime }=m_0/m.$$

(63)

Thus, ${\widehat{P}}_0^{\prime }$ is simply the proportion $m_0/m$ of negative droplets among the m droplets accepted by the droplet reader. Values of both $m_0$ and m are automatically reported by the software that accompanies the ddPCR droplet reader.

Generally speaking, the two best methods for estimating a 95% confidence interval for a binomial success probability like $P_0^{\prime }$, in terms of coverage probability and expected width, are the Wilson method and the Agresti-Coull method [42]. Simulation studies have shown that unless the number of trials is very large ($\gg 100$), the traditional Wald (asymptotic) method typically yields confidence intervals for a binomial probability that correspond to a true confidence level that is much lower than the nominal 95% (meaning that the actual probability that the true value of $P_0$ lies within the interval is much less than the nominal 95%) and hence are much narrower than they should be [42,43]. But because the total number m of accepted droplets is so large in ddPCR applications (typically, $m>$ 10,000), Wald confidence intervals tend to be very similar to the Wilson and Agresti-Coull intervals (see the examples below in Section 3.2.4).

Substituting maximum-likelihood estimator ${\widehat{P}}_0^{\prime }$ for $P_0^{\prime }$, Eqs. (57) and (58) become

$$\begin{array}{lll}\hfill C_{\mathrm{S}}+\mathrm{AM}& =\hfill & -{log}_e(m_0/m)/\overline{v}\end{array}$$

(64)

$$\begin{array}{lll}\hfill C_{\mathrm{sample}}& =\hfill & -{log}_e(m_0/m)F/\overline{v}.\end{array}$$

(65)

These formulas allow us to estimate the unknown concentration of TSC in the S+AM or in a sample, once we have determined the mean droplet volume $\overline{v}$ and, in the latter case, the dilution correction factor F. Note that when applying these formulas, the volume unit of $\overline{v}$ determines the volume unit of $C_{\mathrm{S}}+\mathrm{AM}$ and $C_{\mathrm{sample}}$. Thus, if we wish to express these concentrations in GC/$\mu$L, then $\overline{v}$ must be expressed in $\mu$L instead of nL (note: estimates of $\overline{v}$ usually are reported in nL).

Eqs. (64) and (65) provide point estimates of TSC concentration but give no estimate of their level of uncertainty. As we noted regarding qPCR estimates of TSC copy number per reaction, there are two main sources of this uncertainty: one due to error in the estimate of $\overline{v}$ used, and the other due to statistical error in estimating $P_0^{\prime }$ by $m_0/m$.

Regarding the first source, let $\tilde{\overline{v}}$ and $\overline{v}$ denote the estimated and true values of the mean droplet volume, and suppose $\tilde{\overline{v}}$ contains proportional error ${\eta }_{\mathrm{drop}}\overline{v}$ so that

$$\tilde{\overline{v}}=(1+{\eta }_{\mathrm{drop}})\overline{v}.$$

(66)

Then the estimated sample TSC concentration ${\tilde{C}}_{\mathrm{sample}}$ based on $\tilde{\overline{v}}$ is related to the concentration $C_{\mathrm{sample}}$ based on the correct value $\overline{v}$ as

$${\tilde{C}}_{\mathrm{sample}}={\displaystyle \frac{-{log}_e(m_0/m)F}{\tilde{\overline{v}}}}={\displaystyle \frac{-{log}_e(m_0/m)F}{(1+{\eta }_{\mathrm{drop}})\overline{v}}}={\displaystyle \frac{C_{\mathrm{sample}}}{1+{\eta }_{\mathrm{drop}}}}.$$

(67)

Thus, the proportional error in the estimated TSC concentration in a sample is the inverse of the proportional error in the estimate mean droplet volume. That is,

$${\displaystyle \frac{{\tilde{C}}_{\mathrm{sample}}}{C_{\mathrm{sample}}}}={\displaystyle \frac{1}{1+{\eta }_{\mathrm{drop}}}}.$$

(68)

(The same relationship holds for the ratio of ${\tilde{C}}_{\mathrm{S}}+\mathrm{AM}$ to $C_{\mathrm{S}}+\mathrm{AM}$.) This relationship indicates that overestimates of the mean droplet volume encoded in droplet generator software, as have been well documented in the past, will result in systematic underestimates of sample concentrations. For example, if the original nominal value of 1.00 nL were used fo the mean droplet volume and the actual value is 0.80 nL, then the estimated sample concentration ${\tilde{C}}_{\mathrm{sample}}$ would be about 83% of the correct value, or 17% too low. If an estimate of 0.85 nL were used instead 0.80 nL, ${\tilde{C}}_{\mathrm{sample}}$ would be about 95% of the correct value, or 5% too low. Clearly, it is very important to have an accurate estimate of $\overline{v}$.

The other source or error in ${\tilde{C}}_{\mathrm{sample}}$ can be addressed by obtaining estimates of the lower and upper 95% confidence limits for $C_{\mathrm{S}}+\mathrm{AM}$ or $C_{\mathrm{sample}}$. One way to do this is by substituting the upper and lower 95% confidence limits (respectively) for $P_0^{\prime }$ in Eqs. (64) or (65) (another way is by using bootstrapping [44]). As noted above regarding the confidence limits for $P_0^{\prime }$, the total number m of scanned droplets is large enough so the Wilson, Agresti-Coull, and Wald confidence intervals typically will be very similar. Alternatively, we can estimate confidence limits directly for $C_{\mathrm{S}}+\mathrm{AM}$ or $C_{\mathrm{sample}}$ by bootstrapping; e.g., using the percentile or ${\mathrm{BC}}_a$ algorithm of Efron and Tibshirani [44]. Because m is so large, the confidence limits estimated by these algorithms typically will be very similar to each other, as well as to those based on the confidence intervals for $P_0^{\prime }$ just mentioned (see the examples in Section 3.2.4 below).

These confidence limits do not account for variation in droplet volume, which in any case cannot be quantified unless one can accurately measure the volumes of individual droplets. Currently, the only practical method we know of for estimating individual droplet volumes is the optical microscopy imaging method [34,35,36]. Estimates of the mean and standard deviation of droplet volume obtained in this way can be used to assess the plausibility of treating the sample mean droplet volume $\overline{v}$ as a constant, as well as the plausibility of the related assumption that, to a good approximation, $v_j=\overline{v}$ and hence ${\pi }_j=\overline{\pi }$ for all droplets j that was required in deriving the simple standard formula for TSC concentration C in Section 3.2.1.

Toward this end, let us assume the droplets generated from the S+AM in a particular plate well have true mean volume ${\mu }_v$, and let ${\overline{v}}_m$ denote the sample mean calculated for a set of m droplets whose individual volumes were estimated with optical microscopy imaging. By the strong law of large numbers, ${\overline{v}}_m\to {\mu }_v$ with probability one as $m\to \infty$. But for any finite m, there will be error in estimating ${\mu }_v$ by ${\overline{v}}_m$. Denoting the relative error by $|{\overline{v}}_m-{\mu }_v|/{\mu }_v$, we ask: How likely is it that relative error $|{\overline{v}}_m-{\mu }_v|/{\mu }_v>\delta$, where $\delta$ is a small positive number ($0<\delta <1$). By Chebychev’s inequality,

$$P\left\{|{\overline{v}}_m-{\mu }_v|/{\mu }_v>\delta \right\}\le {\displaystyle \frac{1}{m}}{\left({\displaystyle \frac{\mathrm{CV}}{\delta }}\right)}^2$$

(69)

[e.g., 45, p. 30], where $P\left\{x\right\}$ denotes the probability of event x and CV is the coefficient of variation of droplet volume. Dagata et al. [35] report CV estimates ranging from 0.0042 to 0.0096 for two different DNA sources and two different master mixes. For a human genomic DNA source and a master mix containing dUTP (deoxyuridine triphosphate), they report a CV of 0.0046. Rounding this value up to 0.005, and choosing $\delta =0.1$ and $m=10,000$, we find that

$$P\left\{|{\overline{v}}_m-{\mu }_v|/{\mu }_v>0.1\right\}\le 2.5\times {10}^{-7},$$

(70)

and for a single droplet ($m=1$),

$$P\left\{|v_j-{\mu }_v|/{\mu }_v>0.1\right\}\le 2.5\times {10}^{-3}.$$

(71)

Thus, the coefficient of variation of droplet size is sufficiently small, and the number m of accepted droplets is sufficiently large, so it is reasonable to expect deviations of $v_j$ from ${\mu }_v$ that are more than 1/10-th as large as ${\mu }_v$ to be rare, and deviations of ${\overline{v}}_m$ from ${\mu }_v$ that are more than 1/10-th as large as ${\mu }_v$ to be extremely rare.

An important property of Eq. (64) and (65) is that estimates $C_{\mathrm{S}}+\mathrm{AM}$ and $C_{\mathrm{sample}}$ of TSC concentration depend only indirectly and weakly on the efficiencies of amplification and fluorescence. Because ddPCR is based on an endpoint reaction, amplification and fluorescence only need to be sufficient to eventually produce detectable above-threshold fluorescence in droplets that initially contained one or more TSC. Therefore, in marked contrast to qPCR, ddPCR is relatively insensitive to sample constituents that interfere with amplification or fluorescence, as long as they do not fully block these processes. This property has been confirmed empirically (e.g., [46]).

3.2.3. Estimating the Mean Droplet Volume

Methods for estimating mean droplet volume $\overline{v}$ by optical microscopy imaging are detailed by several authors [34,35,36]. These methods are relatively straightforward and, aside from standard ddPCR instrumentation and supplies, require only a high-quality research microscope with a digital imaging system. Košir et al. [36] found that mean droplet volume showed statistically significant and meaningful differences for different droplet generator models and different master mixes, and they therefore recommended that “to accurately determine copy number, the droplet volume should be measured for all possible droplet generators and master mixes used in a laboratory.” In our experience, however, this admonition is routinely ignored by ddPCR users in the biological and environmental sciences.

A possible alternative method for determining the appropriate mean droplet volume for the specific combination of droplet generator and master mix used in a particular lab is to employ a CRM as a calibration standard. A CRM is simply a sample whose TSC concentration has been estimated by some other method (e.g., chamber-based digital PCR) that is trusted sufficiently so it is thought reasonable to consider the estimated concentration known. Once $m_0$ and m have been determined for a CRM, the only unknown in Eq. (64) and (65) is $\overline{v}$. The appropriate value of $\overline{v}$ for a particular droplet generator and master mix can then be estimated as the value that causes the concentration calculated with Eq. (42) to match the known TSC concentration of a CRM. Rearranging Eq. (65), this value is given by

$$\overline{v}=-{log}_e(m_0/m)F/{\tilde{C}}_{\mathrm{sample}},$$

(72)

where ${\tilde{C}}_{\mathrm{sample}}$ is the estimated TSC concentration in the CRM (expressed in the same volume unit as $\overline{v}$), F is the dilution correction factor, and $m_0$ and m are the number of negative droplets and total number of droplets reported when the CRM is analyzed. Aside from ddPCR instrumentation and supplies, this method requires only the purchase of a CRM. It depends critically on accuracy of the certified concentration, which is also true of standard curves in qPCR analysis. A variety of CRMs are currently available from, for example, the Joint Research Centre’s Institute for Reference Materials and Measurements (JRC-IRMM) and the American Oil Chemists’ Society (AOCS). This method is much simpler and much less time-consuming than the optical microscopy method, but its validity has not yet been assessed.

Once the mean droplet volume for a particular combination of droplet generator and master mix has been determined, TSC concentrations for samples processed with the same droplet generator and master mix can be calculated either by using Eq. (42), with the value of $\overline{v}$ that was determined and the values of $m_0$ and m reported by the instrument software, or by multiplying the concentration reported by the instrument software by the ratio ${\overline{v}}_{\mathrm{soft}}/{\overline{v}}_{\mathrm{meas}}$, where ${\overline{v}}_{\mathrm{soft}}$ is the value of $\overline{v}$ used by the instrument software and ${\overline{v}}_{\mathrm{meas}}$ is the value determined by optical microscopy or calibration with a CRM. If the instrument software does not report the assumed mean droplet volume ${\overline{v}}_{\mathrm{soft}}$, it can be determined from the values of $m_0$, m, and C reported by the instrument software for one or more samples; if these values are inserted on the right side of Eq. (72), the calculated mean droplet volume will be ${\overline{v}}_{\mathrm{soft}}$ (possibly with a small roundoff error).

The effect of error in estimated CRM concentration ${\tilde{C}}_{\mathrm{CRM}}$ is to create error in the estimate of mean droplet volume, which in turn creates error in concentration estimates for field samples. To assess the magnitude of this error, suppose the estimate of CRM concentration ${\tilde{C}}_{\mathrm{sample}}$ (i.e., the certified concentration) contains proportional error $(1+{\eta }_{\mathrm{CRM}})C_{\mathrm{CRM}}$, where $C_{\mathrm{CRM}}$ is the correct CRM concentration. The resulting estimate $\tilde{\overline{v}}$ of mean droplet volume is then

$$\tilde{\overline{v}}={\displaystyle \frac{-log(m_0/m)}{(1+{\eta }_{\mathrm{CRM}})C_{\mathrm{CRM}}}}={\displaystyle \frac{\overline{v}}{1+{\eta }_{\mathrm{CRM}}}}.$$

(73)

In turn, the resulting estimate of sample TSC concentration is

$${\tilde{C}}_{\mathrm{sample}}={\displaystyle \frac{-{log}_e(m_0/m)F}{\tilde{\overline{v}}}}=(1+{\eta }_{\mathrm{CRM}})C_{\mathrm{sample}}.$$

(74)

It follows that

$${\displaystyle \frac{{\tilde{C}}_{\mathrm{sample}}}{C_{\mathrm{sample}}}}=1+{\eta }_{\mathrm{CRM}},$$

(75)

As in the case of qPCR estimates of TSC copy numbers per reaction, the proportional error in calibration standards is passed on unaltered to estimates for field samples, underscoring the importance of minimizing error in estimating the TSC concentration in CRMs. The alternative method of estimating $\overline{v}$ by optical microscopy imaging avoids the issue of CRM concentration accuracy but is much more time consuming and also involves meaurement error.

3.2.4. Example: Droplet Classification and Estimating the TSC Concentration

Figure 8 shows a plot of fluorescence amplitude (vertical axis) versus event or accepted droplet number (horizontal axis) created by Bio-Rad’s ${\mathrm{QuantaSoft}}^{\mathrm{TM}}$ software. In this example, three replicates of a positive control for the human-associated Bacteroidales genetic marker (target sequence) HF183 were analyzed with ddPCR. The horizontal magenta line at amplitude 429 on the vertical axis is the manually-set threshold fluorescence level used to classify droplets as positive (fluorescence amplitude greater than the threshold) or negative (fluorescence amplitude not greater than the threshold). Blue and gray dots are droplets classified as positive and negative, respectively. Sample well labels are listed at the top of the figure; classified droplets for different wells are separated by vertical dashed lines.

The numbers of negative droplets ($m_0$) and all droplets (m) for each of the three wells shown in Figure 8 are shown in Table 4. Also shown are the maximum-likelihood estimates of binomial probability $P_0^{\prime }$ calculated with Eq. (63), the estimated concentration C of the HF183 marker in S+AM calculated with Eq. (64), three types of 95% confidence intervals for $P_0^{\prime }$ (computed with function binom.confint() in R package binom [47]), and six types of 95% confidence intervals for C (bootstrap confidence intervals were computed with function boot.ci() in R package boot [48] and function bcanon() in R package bootstrap [49]). Note that with one exception, the upper and lower 95% confidence limits for C are very similar for all types of confidence interval. The one exception is that the lower 95% confidence limits reported by Bio-Rad’s QuantaSoft software are consistently greater than the lower limits for the other types of confidence interval by more than 1.0  GC/$\mu$L (the QuantSoft limits are based on a Poisson approximation to the actual binomial distribution, but the software documentation does not state which of the many methods for estimating Poisson confidence intervals is used). The mean droplet volume supplied by the instrument manufacturer ($\overline{v}$ = 0.85 nL = $8.5\times {10}^{-4}$ $\mu$L) was used in all calculations because no estimate was available for the specific combination of droplet generator and master mix used.

4. Discussion

Our review and assessment of the theoretical underpinnings of qPCR and ddPCR methods for estimating target sequence copy numbers and concentrations in water samples support the conclusion that both methods are sound. However, careful attention must be paid to specific details to ensure that the estimates obtained are accurate and precise. These key details differ between qPCR and ddPCR.

In the case of qPCR, it is particularly important to reduce any significant PCR inhibition by sample constituents, since qPCR copy number estimates depend strongly on the rate of amplification (amplicons produced per DAE cycle) and therefore are very sensitive to amplification efficiency. It is also important to ensure proper choice and fitting of the regression model for the standard curve, including assessment of potential outliers and potential nonlinearity. If one wishes to estimate confidence intervals for the slope or intercept parameters (e.g., for quality assurance purposes) or estimate confidence intervals for copy-number estimates produced by inverse prediction, then the statistical assumptions of the regression model must be carefully assessed. If any of these assumptions are not tenable, an alternative model must be chosen and assessed. This model selection process is nontrivial for three main reasons:

• Over the range of standards typically required for analysis of water samples, the variance of measured $c_{\mathrm{T}}$ values often differs markedly for different standards, meaning that the $c_{\mathrm{T}}$ values are heteroskedastic. In such cases, the variance homogeneity assumption of classical OLS regression is not tenable.
• One way to address heteroskedasticity is by employing WLS regression. However, any choice of weights must be carefully assessed to ensure it succeeds in homogenizing the variance in $c_{\mathrm{T}}$ residuals.
• Even if weights can be found so that WLS regression successfully homogenizes the variance of residuals, the resulting intercept and slope parameters may exaggerate errors in predicted sample copy numbers at low concentrations, due to heavier weighting of residuals for high concentrations, where $c_{\mathrm{T}}$ values typically are less variable.

In the case of ddPCR, an important but commonly ignored issue is that the mean droplet volume—to which concentration estimates are quite sensitive—is known to differ meaningfully for different combinations of droplet generator and master mix. The fact that mean droplet volume is not constant across different instruments and assay mixes implies that no universal estimate of this critical parameter (e.g., a single value coded into instrument software by the instrument manufacturer) should be trusted. Instead, it is advisable to determine the appropriate mean droplet volume for the specific combination of droplet generator and master mix being used for a particular study, either by optical microscopy or by calibration with a CRM.

Both qPCR and ddPCR are potentially appropriate for many of the same types of analysis in studies of aquatic systems and may therefore be viewed as alternatives that one must choose between. With this choice in mind, we now summarize what we see as the main advantages and disadvantages of qPCR and ddPCR. Our comparison of methods is structured around six groups of properties: cost, sample turnaround time, calibration, inhibition, limits of quantification, and simplicity (Table 5). Basu [9] provides an alternative comparison of these methods (which we recommend to the reader) that considers some of the same properties but several others, as well. For the properties that both comparisons consider, our views regarding the advantages and disadvantages of qPCR and ddPCR largely agree with those of Basu [9].

Cost. The instrument and per-sample costs for ddPCR are substantially higher than those for qPCR at the time of this writing, but we expect both types of cost to become more similar to those for qPCR in the future.

Sample turnaround time. The difference in time required to prepare and analyze samples is mainly important in applications where decisions based on the results must be made as soon as possible after sample collection, as is often the case when monitoring fecal indicator bacteria at recreational beaches to assess compliance with a Beach Action Value for protecting human health. In our experience with water quality monitoring and microbial source tracking using the ${\mathrm{StepOnePlus}}^{\mathrm{TM}}$ real-time qPCR system (Applied ${\mathrm{Biosystems}}^{\circledR }$) and the ${\mathrm{QX200}}^{\mathrm{TM}}$ ddPCR system (${\mathrm{Bio-Rad}}^{\circledR }$), the turnaround time is roughly two hours longer for ddPCR (5–6 h) than for qPCR (3–4 h). This difference is due mainly to the extra time required by ddPCR to read and count droplets from each well of a 96-well plate, which is done sequentially rather than in parallel.

Calibration. Contra an opinion that is commonly expressed in the literature, both qPCR and ddPCR require some form of calibration to achieve acceptable accuracy. For qPCR, this involves fitting a standard curve to a set of calibration data for each set of analyses performed, as well as periodic purchasing and running of calibration plates provided by the instrument manufacturer. For ddPCR, a calibration procedure to determine the mean droplet volume for the specific combination of droplet generator and master mix used for a particular set of analyses is strongly recommended [34,35,36]; any value of this crucial parameter that was not estimated for the specific combination of droplet generator and master mix employed in a particular study is likely to differ meaningfully from the actual mean droplet volume, and if it does, the result will be systematic and meaningful errors in estimates of target sequence copy numbers or concentrations in samples.

Inhibition. ddPCR is inherently less sensitive to PCR inhibition than is qPCR, mainly because it is an endpoint PCR method and therefore depends only weakly on amplification efficiency [50,51,52]. This is an important advantage of ddPCR when significant inhibition by sample constituents occurs, which sometimes happens in studies of aquatic systems. In some cases, inhibition in qPCR analyses can be adequately reduced by sample dilution, sample cleaning, and/or use of a master mix such as TaqMan™ Environmental Master Mix 2.0 that is specifically formulated for this purpose [52].

Limits of quantification. Several studies have found that the lower limit of quantification is lower for ddPCR than for qPCR [51,53,54], two have found that the upper limit of quantification without sample dilution is higher for qPCR [52,54], and one study notes that the dynamic range (difference between the upper and lower limits of quantification) is substantially wider for qPCR [9]. In our opinion, however, these patterns require further assessment with a greater variety of sample types, including field samples containing PCR inhibitors. It is also important to bear in mind that one can relatively easily and quickly increase the effective upper limit of quantification by diluting samples, while decreasing the lower limit of quantification by concentrating samples is more time-consuming and may not be feasible.

Simplicity. It is our experience, based on working with many students, that qPCR is somewhat easier to perform properly than is ddPCR. On the other hand, the methods of data analysis and their theoretical underpinnings are much simpler for ddPCR than for qPCR, with the main difference from a data analyst’s perspective being that the ticklish statistical issues related to proper choice and fitting of a regression model for the standard curve are avoided when ddPCR is used.

Overall, the facts that ddPCR does not require a standard curve and therefore avoids the statistical issues regarding proper fitting and assessment of the curve, is inherently less sensitive to PCR inhibitors than is qPCR, and has lower limits of quantification all favor the choice of ddPCR. This is especially true in basic research studies of aquatic systems where a significant proportion of the field samples contain low concentrations of the target sequence and/or exhibit strong PCR inhibition that cannot be adequately reduced, and where high levels of data censoring are not acceptable. On the other hand, the significantly lower costs associated with qPCR, the shorter turnaround time for obtaining results, the higher upper limit of quantification, and the wider dynamic range all favor qPCR when funding levels are low, the concentrations that are most important to quantify accurately are well within the method’s range of quantification (as is often true in monitoring studies), and strong inhibition is not present or can be adequately reduced, especially when a short turnaround time for obtaining results is required.

In our opinion, qPCR often will be preferable in laboratory studies where investigators have significant control over the range of target sequence concentrations in samples, the concentrations are sufficiently high so the set of appropriate standards does not include concentrations that are low enough to exhibit exaggerated variability in $c_{\mathrm{T}}$, and no significant inhibition is present. The choice of method in field studies is more likely to be application-specific. The advice of Chandler [55] still applies today: “To realize the full potential of the PCR in environmental situations, ... one should be aware of its limitations at levels of chemical and genetic complexity not normally encountered in traditional molecular biology laboratories.” Several of those limitations differ between qPCR and ddPCR, and they are the ones that will determine which method makes more sense for a particular environmental application.

5. Conclusions

• The theoretical basis of mathematical and statistical methods commonly used for estimating target sequence copy numbers and concentrations with qPCR and ddPCR is sound.
• The reliance of qPCR on a standard curve creates both complications and uncertainties in fitting and assessing the standard curve, because the calibration data typically are heteroskedastic.
• Compared to ddPCR, the method for estimating copy numbers and concentrations with qPCR is more sensitive to sample properties that interfere with fluorescence intensity or reduce amplification efficiency, making the use of effective methods to reduce interference particularly important.
• Estimating copy numbers and concentrations with ddPCR does not rely on a standard curve and therefore avoids statistical complications and uncertainties regarding the proper fitting and assessment of standard curves when the calibration data are heteroskedastic.
• Accuracy of ddPCR copy number and concentration estimates is sensitive to the mean droplet volume, which differs meaningfully for different combinations of droplet generator and master mix. Therefore, the mean droplet volume should be determined empirically for the particular combination of droplet generator and master mix used in a given analysis instead of relying on a rough universal estimate (e.g., one coded into software supplied by the instrument manufacturer).