Estimation of DNA copy number in a given biological sample is an extremely important problem in genomics. This problem is especially challenging when the number of the DNA strands is minuscule, which is often the case in applications such as pathogen and genetic mutation detection. A recently developed technique, real-time polymerase chain reaction (PCR), amplifies the number of initial target molecules by replicating them through a series of thermal cycles. Ideally, the number of target molecules doubles at the end of each cycle. However, in practice, due to biochemical noise the efficiency of the PCR reaction, defined as the fraction of target molecules which are successfully copied during a cycle, is always less than 1. In this paper, we formulate the problem of joint maximum-likelihood estimation of the PCR efficiency and the initial DNA copy number. As indicated by simulation studies, the performance of the proposed estimator is superior with respect to competing statistical approaches. Moreover, we compute the Cramer-Rao lower bound on the mean-square estimation error

Hassibi, A.

Hassibi, B.

Vikalo, H.

English

Caltech Authors - Main

ML ESTIMATION OF DNA INITIAL COPY NUMBER IN POLYMERASE CHAIN
REACTION (PCR) PROCESSES
H. Vikalo∗, B. Hassibi∗, and A. Hassibi‡
∗Electrical Engineering Department, California Institute of Technology, Pasadena, CA
‡Electrical and Computer Engineering Department, University of Texas, Austin, TX
ABSTRACT
Estimation of DNA copy number in a given biological
sample is an extremely important problem in genomics.
This problem is especially challenging when the num-
ber of the DNA strands is minuscule, which is often the
case in applications such as pathogen and genetic muta-
tion detection. A recently developed technique, real-time
polymerase chain reaction (PCR), ampli?es the number
of initial target molecules by replicating them through
a series of thermal cycles. Ideally, the number of tar-
get molecules doubles at the end of each cycle. How-
ever, in practice, due to biochemical noise the ef?ciency
of the PCR reaction, de?ned as the fraction of target mole-
cules which are successfully copied during a cycle, is al-
ways less than 1. In this paper, we formulate the problem
of joint maximum-likelihood estimation of the PCR ef?-
ciency and the initial DNA copy number. As indicated by
simulation studies, the performance of the proposed esti-
mator is superior with respect to competing statistical ap-
proaches. Moreover, we compute the Cramer-Rao lower
bound on the mean-square estimation error.
Index terms – biomedical signal processing, estimation,
error analysis
1. INTRODUCTION
Ampli?cation and quanti?cation of DNA strands with poly-
merase chain reaction (PCR) process [1] is an essential
part of many biotechnological procedures. Applications
of PCR [2] include genotyping, detection of infectious
and hereditary diseases, and genetic ?ngertyping, to name
a few. Typically, a given biological sample contains only
a small amount of a target DNA. PCR ampli?es the tar-
get in vitro, relying on an enzymatic replication process in
each of its temperature-regulated cycles (typically, 30−40
of them). A PCR cycle, illustrated in Figure 1, consists
of three distinct phases: denaturing, annealing, and ex-
tension. During denaturing, the sample is heated (typi-
cally, to temperatures above 90◦C) to break the hydro-
gen bonds between strands of the target DNA fragments,
creating twice as many single-stranded fragments. Each
This work was supported in part by the National Science Founda-
tion under grant no. CCR-0133818, by the Of?ce of Naval Research
under grant no. N00014-02-1-0578, and by Caltech’s Lee Center for
Advanced Networking.
of the single-stranded DNA fragments serves as a tem-
plate during the second phase, annealing, when the sam-
ple is cooled down to temperatures typically between 40−
72◦C. At such temperatures, primers – short sequences
of nucleotides, designed to be exact complements to spe-
ci?c regions on the templates – hybridize to the templates.
There are two types of primers, one for each of the two
types of templates. The primers serve as initiation sites
for a DNA polymerase enzyme activated in the last phase
of a cycle, extension. The sample is heated to 72◦C max-
imizing the rate of extension while ensuring that the par-
tially extended primers remain attached to the templates.
Ideally, at the end of the extension phase, there are twice
as many double-stranded target DNA strands as there were
at the beginning of the cycle. This implies an exponential
growth of the number of the target DNA. However, prac-
tical issues affect the replication process adversely and
the ef?ciency of PCR – de?ned as the probability of gen-
erating a replica of each template molecule – is smaller
than 1. The random nature of the underlying biochemical
process leads to variations in the PCR yield. Moreover,
the creation of non-speci?c byproducts in the replication
process further diminishes the purity of the PCR product.
The probabilistic nature of the replication process is ad-
dressed in [6]-[9], where various stochastic models have
been proposed. In [10], the mutations-related effects that
plague the ef?ciency of PCR have been studied.
Fig. 1. Typical PCR cycle
Quanti?cation of the created amplicons (DNA mole-
cules obtained by the replication of the initial DNA strand)
is based on measuring the light intensity originating from
the ?uorescent reporter molecules incorporated into the
I ­ 4171­4244­0728­1/07/$20.00 ©2007 IEEE ICASSP 2007
double-stranded DNA (dsDNA). One such reporter is SYBR
Green I, a dye that binds to dsDNA after which its ?uo-
rescence activity increases signi?cantly. Other reporters
include hybridization and TaqMan probes (see, e.g., [3]).
In real-time PCR, the ?uorescent signal is measured
at the end of each cycle. The measured light intensities
comprise a reaction pro?le, usually plotted against the
number of cycles. A typical reaction pro?le has three dis-
tinct phases: a background phase, an exponential growth
phase, and a saturation phase. During the ?rst phase,
the background noise originating from unbound probes
dominates the useful signal emanating from the probes
attached to the templates. Although the ?uorescent level
of the unbound probes is much lower than the ?uorescent
level of the probes bounded to the double-stranded target
DNA, the former signi?cantly outnumber the latter dur-
ing the ?rst 15− 20 cycles. The second phase starts when
the signal from the PCR products rises suf?ciently above
the background noise. Typically, measurements collected
during the exponential growth phase (also referred to as
the log phase) are the only ones used to infer information
about the original number of the DNA targets in the bi-
ological sample. The reason for imposing such a restric-
tion is that the ef?ciency of PCR can be assumed con-
stant in the ?rst two phases, which makes the estimation
tractable. In the third phase of PCR, however, the ef?-
ciency decreases rapidly as the reaction enters a plateau.
The ultimate goal of PCR is the estimation of the ini-
tial number of target molecules. In practice, this is com-
monly done by comparing a PCR reaction pro?le with the
reaction pro?le of a so-called standard, where the latter
are recorded for several initial concentrations of a target
which has the same ef?ciency as the DNA target of inter-
est. In recent work [7], [8], the reaction pro?le is used
directly (i.e., without use of a standard) to estimate the ef-
?ciency, which in turn is then used to ?nd an estimate of
the initial number of the DNA target molecules.
In this paper, we ?nd the joint maximum-likelihood
estimate of the PCR ef?ciency and the number of initial
target molecules. Furthermore, we ?nd the Cramer-Rao
lower bound on the minimum mean-square error of the
estimated parameters and illustrate by simulations that the
proposed estimator can achieve it.
2. THE MODEL AND JOINT ML ESTIMATION
Let x0 denote the initial copy number of target DNA mole-
cules which we want to estimate. We assume that the ef-
?ciency of replication during both the background phase
and the exponential phase is constant, and denote it by p.
(During the saturation phase, the ef?ciency drops as the
reaction approaches a plateau. For the sake of simplic-
ity of the estimation procedure, we use only the measure-
ments taken at the end of the cycles wherein the ef?ciency
is constant.) Furthermore, denote the number of target
molecules at the end of the nth cycle by xn, and note that
xn = xn−1 + an,
where an is the number of amplicons that have been cre-
ated in the nth cycle. Since the probability that each of
the xn−1 available amplicons extends in the nth cycle is
p, it is easy to see that an is a binomial random variable
with mean pxn−1 and variance p(1 − p)xn−1. We may
therefore write
xn = (1 + p)xn−1 + x˜n, (1)
where x˜n is a random variable with zero mean and vari-
ance p(1 − p)xn−1. Recursion (1) describes a branching
process, often used to model replication in biological sys-
tems [5]. It is not too dif?cult to show (see, e.g., [6]) that
the mean of xn in (1) is given by
E{xn} = (1 + p)nx0. (2)
Furthermore, its variance can be found as ([6])
σ2n =
1− p
1 + p
[
(1 + p)2n − (1 + p)n]x0. (3)
Imperfect instrumentation and other biochemistry in-
dependent sources create a noise which corrupts the mea-
surements of xn. We assume that the noise is additive
Gaussian N (0, σ2w), and denote it by wn. Hence, the
quantity measured is given by
zn = xn + wn.
Let us denote the number of temperature cycles in the
background phase of PCR by k. Therefore, the ?rst mea-
surement taken beyond the background noise level is zk+1.
Furthermore, denote the number of temperature cycles in
the exponential phase by l. Hence, the last measurement
taken before the ef?ciency starts to rapidly deteriorate is
zk+l. Introduce a new variable, y, de?ned as
y =
⎡
⎢⎢⎢⎢⎢⎣
zk+1−(1+p)k+1x0
σk+1
zk+2−(1+p)k+2x0
σk+2
...
zk+l−(1+p)k+lx0
σk+l
⎤
⎥⎥⎥⎥⎥⎦
.
Note that y is zero-mean. Finding the exact probability
density function (pdf) of y appears to be dif?cult. On the
other hand, we can express y as a sum of x0 identically
distributed random variables,
y = y1 + y2 + · · ·+ yx0 ,
where each summand is given by
yi =
⎡
⎢⎢⎢⎢⎢⎣
xi,k+1+wk+1/x0−(1+p)k+1
σk+1
xi,k+2+wk+2/x0−(1+p)k+2
σk+2
...
xi,k+l+wk+l/x0−(1+p)k+l
σk+l
⎤
⎥⎥⎥⎥⎥⎦
,
I ­ 418
where xi,j denotes the number of the amplicons at the end
of the jth cycle which have originated from the ith target
molecule in the original sample. Moreover, y1,y2, . . . ,yx0
are independent since they are comprised of the measure-
ments taken in the ?rst two phases of PCR, during which
targets do not inhibit each other’s replication. Therefore,
y can be represented as a sum of a large number of zero-
mean independent, identically distributed (iid) random vari-
ables. We thus invoke the central limit theorem to argue
that the distribution of y may be approximated by a multi-
variate Gaussian distribution.
Note that the (i, j)-entry of the l× l covariance matrix
of y, R, is given by
R(i, j) = E{yiyj}
=
E{zk+izk+j} − (1 + p)2k+i+jx20
σk+iσk+j
, (4)
where yi and yj denote the ith and the jth component of
y, respectively. We can show that
R(i, j) = (1 + p)j−i
σk+i
σk+j
+
σ2w
σk+iσk+j
δi−j ,
where
δi−j =
{
1, if i = j,
0, otherwise,
Now that we computed the covariance matrix R, the
probability density function of y can be approximated by
the multi-variate Gaussian distribution
fy(y) =
1
(2π)l/2(detR)1/2
e−
1
2y
T R−1y. (5)
Note that fy(y) depends on x0 and p through both y and
R.
2.1. Optimal estimation of x0 and p
The joint maximum-likelihood estimate of x0 and p can
be found by solving the maximization problem
max
x0,p
fy(y),
or, equivalently, by solving the minimization
min
x0,p
{yT R−1y + log detR}. (6)
On the other hand, the traditional approach to the es-
timation of the initial population in a branching process
(see, e.g., [11]) ?rst focuses on ?nding the maximum-
likelihood estimator of p,
pˆ =
zk+1 + · · ·+ zk+l
zk + · · ·+ zk+l−1 − 1. (7)
Then, the above estimate pˆ is used to estimate x0 as
xˆ0 =
zk+l
(1 + pˆ)k+l
. (8)
Note that for the reliability of the estimate pˆ in (7), we
only use measurements taken in the exponential phase of
the PCR.
The objective function of the optimization (6) is not
convex. To solve it, one can use, e.g., a gradient search
initialized by xˆ0 and pˆ obtained from (8) and (7), respec-
tively.
3. LIMITS OF PERFORMANCE OF PCR: THE
CRAMER-RAO BOUND
The minimum mean-square error of any estimation pro-
cedure is lower bounded by the Cramer-Rao bound [12].
We will derive the Cramer-Rao lower bound (CRLB) and
use it to quantify the limits of achievable performance of
DNA estimation in PCR.
Collect the parameters that need to be estimated into
a vector, cT = [x0 p]
T . The Fisher information matrix,
F , is given by the negative of the expected value of the
Hessian matrix of log py|c(y), i.e.,
F = −Ey
{∇c∇Tc log py|c(y)} .
Therefore, the entries of the 2× 2 matrix F are given by
Fij = −Ey
{
∂2
∂ci∂cj
log py|c(y)
}
,
where, for compactness of the notation, ci and cj denote
the entries of c (i.e., c1 = x0, c2 = p). Assuming an unbi-
ased estimator, the CRLB on the minimum mean-square
error of estimating x0 is given by
E
{
(xˆ0 − x0)2
} ≥ [F−1]11, (9)
where [F−1]11 denotes the (1, 1)-entry of F−1. Similarly,
the CRLB on the minimum mean-square error of estimat-
ing p is
E
{
(pˆ− p)2} ≥ [F−1]22, (10)
where [F−1]22 denotes the (2, 2)-entry of F−1.
Let us denote L1(c) = log detR and L2(c) = yT R−1y,
so that we can write
L(c) = log(py|c(y)) = − log(2π)− 12L1(c)−
1
2
L2(c).
Therefore, the Fisher information matrix can be written
as
F =
1
2
Ey
{∇c∇Tc L1(c)}+ 12Ey
{∇c∇Tc L2(c)} .
I ­ 419
It is easy to ?nd the components of ∇c∇Tc L1(c),
∂2
∂ci∂cj
L1(c) = Tr
{
−R−1 ∂R
∂cj
R−1
∂R
∂ci
+ R−1
∂2R
∂cj∂ci
}
,
(11)
where Tr {·} denotes the trace operation over its argu-
ment. Finding the components of ∇c∇Tc L2(c) is some-
what more involved. Due to space limitation, we omit the
derivation – although straightforward, the ?nal expres-
sions for Fij is fairly cumbersome to write (for details,
we refer the interested reader to [13]). It suf?ces to say
that one may use a symbolic math manipulation package
(e.g., Mathematica, Maple) to ef?ciently compute Fij for
a given set of parameters (σ2w, k, l, x0, p). This, in fact,
is how we proceed: we use Mathematica to evaluate the
CRLBs in (9) and (10) for any given set of the RT-PCR
experiment parameters.
4. SIMULATION RESULTS AND CONCLUSION
In Figure 2, we compare the mean-square error of the es-
timate of x0 computed by (6) and that of (8), and compare
them with the Cramer-Rao lower bound. The PCR is sim-
ulated as a branching process with x0 = 1000, while the
variance of the noise in the exponential phase is assumed
to be 1/100 of the measured signal intensity.
0.75 0.8 0.85 0.9 0.95
0
50
100
150
200
250
300
350
efficiency p
sq
u
a
re
 r
o
o
t 
o
f 
m
e
a
n
?
sq
u
a
re
 e
rr
o
r
x
0
 from the joint ML estimator
simple estimator x
n+1
/(x
n+1
/x
n
?1)n+1
Cramer?Rao bound
Fig. 2. Comparison of the estimation mean-square errors
E(xˆ0 − x0)2
From Figure 2, we see that the proposed joint ML es-
timator (6) signi?cantly outperforms the estimator (8) for
all considered values of p. Furthermore, the mean-square
error of the joint ML estimator is almost achieving the
Cramer-Rao lower bound. The slight discrepancy could
be caused by the approximation of the distribution of y
by a Gaussian.
It is of interest to extend the results presented here to
the case where the ef?ciency is not constant, but changes
according to a known model. Furthermore, a study of the
PCR in the saturation phase is also of interest.
5. REFERENCES
[1] K. Mullis and F. Faloona, “Speci?c synthesis of
DNA in vitro via a polymerase-catalyzed chain re-
action,” Methods Enzymol 155:335-350, 1987.
[2] M. A. Innis, D. H. Gelfand, and J. J. Sninsky, PCR
Applications: Protocols for Functional Genomics,
Academic Press, 1999.
[3] C.T. Wittwer, M. G. Herrmann, A.A. Moss, and R.P.
Rasmussen, “Continuous ?uorescence monitoring
of rapid cycle DNA ampli?cation,” BioTechniques,
22:130-138, 1997.
[4] P. M. Holland, R. D. Abramson, R. Watson, and
D. H. Gelfand, “ Detection of Speci?c Polymerase
Chain Reaction Product by Utilizing the 5” to 3”
Exonuclease Activity of Thermus aquaticius DNA
polymerase,” PNAS, 88, 7276-80, 1991.
[5] M. Kimmel and D. E. Axelrod, Branching Processes
in Biology, Springer-Verlag, New York, 2002.
[6] G. Stolovitzky and G. Cecchi, “Ef?ciency of DNA
replication in the polymerase chain reaction,” PNAS,
vol. 93, pp. 12947-12952, November 1996.
[7] C. Jacob and J. Peccoud, “Estimation of the parame-
ters of a branching process from migrating binomial
observation,” Adv. in Appl. Prob., 30,948-967, 1998.
[8] N. Lalam, C. Jacob, and P. Jagers, “Modelling
the PCR ampli?cation process by a size-dependent
branching process and estimation of the ef?ciency,”
Adv. in Appl. Prob., 36, 602-615, 2004.
[9] A. Hassibi, H. Kakavand, and T. H. Lee, “A stochas-
tic model and simulation algorithm for polymerase
chain reaction (PCR) systems,” GENSIPS 2004.
[10] D. Wang, C. Zhao, R. Cheng, F. Z. Sun, ”Esti-
mating the mutation rate during error-prone poly-
merase chain reaction,” J. of Computational Biol-
ogy, 7, 143-158, 2000.
[11] J.-P. Dion, “Estimation of the mean and the initial
probabilities of a branching process,” J. of Applied
Probability, 11, 687-694, 1974.
[12] H. Cramer, Mathematical Models of Statistics,
Princeton University Press, Princeton, NJ 1946.
[13] H. Vikalo, B. Hassibi, and A. Hassibi, ”ML
estimation of DNA initial copy number in
polymerase chain reaction (PCR) process,”
http://www.its.caltech.edu/˜ hvikalo/PCRest.pdf.
I ­ 420


ML Estimation of DNA Initial Copy Number in Polymerase Chain Reaction (PCR) Processes

Crossref

A stochastic model and simulation algorithm for polymerase chain reaction (PCR) systems,” GENSIPS

Branching Processes in Biology,

Continuous ƀuorescence monitoring of rapid cycle

Efſciency of DNA replication in the polymerase chain reaction,”

Estimating the mutation rate during error-prone polymerase chain reaction,”

Estimation of the mean and the initial probabilities of a branching process,”

Estimation of the parameters of a branching process from migrating binomial observation,” Adv.

Mathematical Models of Statistics,

ML estimation of DNA initial copy number in polymerase chain reaction (PCR) process,” http://www.its.caltech.edu/˜ hvikalo/PCRest.pdf.

Modelling the PCR ampliſcation process by a size-dependent branching process and estimation of the efſciency,” Adv.

PCR Applications: Protocols for Functional Genomics,

Speciſc synthesis of DNA in vitro via a polymerase-catalyzed chain reaction,” Methods Enzymol

http://authors.library.caltech.edu/9072/1/VIKicassp07.pdf

ML Estimation of DNA Initial Copy Number in Polymerase Chain Reaction (PCR) Processes

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main

Crossref