JP2007226771A - Pcr elbow determination using curvature analysis of double sigmoid - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide systems and methods for determining characteristic transition values such as elbow values in sigmoid or growth-type curves, such as the cycle threshold (Ct) value in PCR amplification curves. <P>SOLUTION: A double sigmoid function with parameters determined by a Levenberg-Marquardt (LM) regression process is used to find an approximation to a curve that fits a PCR dataset. Once the parameters have been determined, the curve can be normalized using one or more of the determined parameters. After normalization, the normalized curve is processed to determine the curvature of the curve at some or all points along the curve, e.g., to produce a dataset or plot representing the curvature v. the cycle number. The cycle number at which the maximum curvature occurs corresponds to the Ct value. The Ct value is then returned and may be displayed or otherwise used for further processing. <P>COPYRIGHT: (C)2007,JPO&INPIT

Description

  The present invention relates generally to systems and methods for processing data representing sigmoid or growth curves, and more particularly to characteristic cycle threshold (Ct) or elbow values of PCR amplification curves, or others The present invention relates to a system and method for determining the elbow value of a growth curve.

  Polymerase chain reaction (PCR) is an extracorporeal method in which a specific nucleic acid sequence is synthesized or amplified by an enzyme. This reaction typically uses two oligonucleotide primers that hybridize to the opposite strand and flanking the template or target DNA sequence to be amplified. The extension of the primer is catalyzed by a thermostable DNA polymerase. Specific DNA fragments accumulate exponentially as a result of a series of repeated cycles involving template denaturation, primer annealing, and polymerase extension of the annealed primer. This process typically uses fluorescent probes or markers to facilitate the detection and quantification of the amplification process.

  FIG. 1 shows a typical real-time PCR curve, where fluorescence intensity values for a typical PCR process are plotted against cycle number. In this case, the formation of the PCR product is monitored at each cycle of the PCR process. Amplification is typically measured in a thermocycler that includes components and devices that measure the fluorescence signal during the amplification reaction. An example of such a thermocycler is Roche Diagnostics LightCycler (Cat. No. 201110468). For example, the amplification product is detected by a fluorescently labeled hybridization probe that emits a fluorescent signal only when bound to the target nucleic acid, or in some cases by a fluorescent dye that binds to double-stranded DNA.

  In a typical PCR curve, identifying the transition point at the end of the baseline region (which is commonly referred to as the elbow value or cycle threshold (Ct)) will understand the characteristics of the PCR amplification process. Very useful to. The Ct value can be used as a measure of the efficiency of the PCR process. For example, a defined signal threshold is typically determined for all reactions to be analyzed, and the number of cycles (Ct) required to reach this threshold is determined for a target nucleic acid and a reference nucleic acid such as a standard or housekeeping gene. . Based on the Ct values obtained for the target nucleic acid and the reference nucleic acid, the absolute or relative copy number of the target molecule can be determined ("Genome Research" by Gibson et al., Vol. 6, pages 995-1001), Bieche. "Cancer Research" (59, 2759-2765, 1999), WO 97/46707, WO 97/46712, WO 97/46714). The elbow value of the region 20 at the end point of the baseline region 15 in FIG. 1 can be located in the region of 30 cycles.

  The elbow value of the PCR curve can be determined using several existing methods. For example, in various current methods, the actual value of the elbow is determined as the value at which the fluorescence reaches a predetermined level called AFL (Arbitrary Fluorescence Value). Other current methods may use the cycle number at which the second derivative of the fluorescence versus cycle number reaches a local maximum. Both of these methods have drawbacks. For example, some methods are very sensitive to outlier (noise) data, and the AFL value method does not work well for data sets with a high baseline. The conventional method of determining the baseline stop point (or baseline endpoint) of the growth curve shown in FIG. 1 cannot work satisfactorily, especially in high titer environments. Furthermore, these algorithms are usually not well defined, are linearly dependent, and many (eg, 50 or more) that are very difficult to optimize if not impossible. ) Parameters.

  Accordingly, it would be desirable to provide a system and method for determining the elbow value of a curve (especially a PCR curve) such as a sigmoid type or growth curve that overcomes the aforementioned and other problems.

  The present invention provides a new and efficient system and method for determining characteristic transition values such as elbow values of sigmoid or growth type curves. In certain implementations, the systems and methods of the present invention are particularly useful for determining the cycle threshold (Ct) of a PCR amplification curve.

  In accordance with the present invention, an approximation to a curve that fits a PCR data set is detected by using a double sigmoid with parameters determined by a Levenberg-Marquardt (LM) regression process. Once the parameter determination is complete, one or more of the determined parameters can be used to normalize the curve. After normalization, the normalized curve is processed to determine the curvature of the curve at some or all points along the curve, for example, generating a data set or plot that represents curvature versus number of cycles. The number of cycles in which the maximum curvature occurs corresponds to the Ct value. This returns a Ct value that can be displayed or used for further processing.

In a first aspect of the invention, a computer-implemented method for determining a point at the end of a baseline region of a growth curve comprising:
-Receiving a data set representing a growth curve, each data set comprising a plurality of data points having a pair of coordinate values;
Determining the parameters of the function by applying a Levenberg-Marquardt (LM) regression process to the double sigmoid function to calculate an approximation of the curve that fits the data set;
Normalizing the curve using the determined parameters to generate a normalized curve;
Determining a maximum curvature point representing the end point of the baseline region of the long curve by processing the normalized curve;
A method comprising the steps of:

In a second aspect of the invention, a computer readable medium is provided that includes code for controlling a processor to determine a point at the end of a baseline region of a growth curve, the code comprising:
Receiving a data set representing a growth curve, the data set comprising a plurality of data points each having a pair of coordinate values;
Determining the parameters of the function by applying a Levenberg-Marquardt (LM) regression process to the double sigmoid function to calculate an approximation of the curve that fits the data set;
Normalizing the curve using the determined parameters to generate a normalized curve;
Determining a point of maximum curvature by processing a normalized curve, the point of maximum curvature representing the end of the baseline region of the growth curve; and
Contains instructions to execute.

In yet another aspect of the invention, a kinetic polymerase chain reaction (PCR) system comprising:
A kinetic PCR analysis module that generates a PCR data set representing a kinetic PCR amplification curve, the data set comprising a plurality of data points each having a pair of coordinate values; and The data set includes a kinetic PCR analysis module that includes data points within a region of interest that includes a cycle threshold (Ct);
-PCR data set
Determining the parameters of the function by applying a Levenberg-Marquardt (LM) regression process to the double sigmoid function to calculate an approximation of the curve that fits the data set;
Normalizing the curve using the determined parameters to generate a normalized curve;
Determining a point of maximum curvature representing the cycle threshold (Ct) of the growth curve by processing the normalized curve;
An intelligence module adapted to process and determine a Ct value;
have.

  Reference to the remaining portions of the specification, including the attached drawings and claims, will realize other features and advantages of the present invention. Further features and advantages of the present invention, as well as the structure and operation of various embodiments of the present invention, are described in detail below in connection with the accompanying drawings. In the accompanying drawings, the same or functionally similar elements are indicated by similar reference numerals.

  The present invention provides a system and method for determining the endpoint of a baseline region of a PCR amplification curve or the transition value of a sigmoid or growth curve, such as an elbow value or Ct value. In one aspect, an approximation to the curve is detected by using a double sigmoid function having parameters determined by a Levenberg-Marquardt (LM) regression process. Once the parameter determination is complete, the curve can be normalized using one or more of the determined parameters. After normalization, the normalized curve is processed to determine the curvature of the curve at some or all points along the curve, for example, to generate a data set or plot that represents curvature versus number of cycles. Is done. The number of cycles in which the maximum curvature occurs corresponds to the Ct value. This returns a Ct value that can be displayed or used for further processing.

  As an example, a growth or amplification curve 10 in the environment of a PCR process is shown in FIG. As illustrated, the curve 10 includes a lag phase region 15 and a log phase region 25. The lag phase region 15 is generally called a baseline or a baseline region. Such a curve 10 includes a transition region 20 of interest that connects the lag phase region and the log phase region. Region 20 is commonly referred to as an elbow or elbow region. The elbow region typically defines the endpoint of the baseline and the transition in the underlying process growth or amplification rate. Identifying specific transition points within region 20 may be useful for analyzing the behavior of the underlying process. In a typical PCR curve, identifying transition points called elbow values or cycle thresholds (Ct) is very useful in understanding the efficiency characteristics of the PCR process.

  Other processes that can provide similar sigmoid or growth curves include bacterial processes, enzymatic processes, and binding processes. For example, in the bacterial growth curve, the transition point of interest is referred to as the time point θ in the lag phase. Other specific processes for generating a data curve that can be analyzed in accordance with the present invention include the SDA (Strand Displacement Amplification) process, the NASBA (Nucleic Acid Sequence-Based Amplification) process, and the TMA (Transcribing Meditation Amplification process). . For examples of SDA and NASBA processes and data curves, see "Homogeneous Real-Time Detection of Single-Nucleotide Solid Strength Disposition by Tang, Sha-Sha, et al." 49 (10), p. 1599), Weusten, Jos. A. M.M. Others "Principles of Quantification of Virtual Loads Using Nucleic Acid Sequence 26 Based A s Accordingly, in the remainder of this specification, examples and embodiments of the present invention will be described in terms of their applicability to PCR curves, but the present invention is based on other process related data curves. It should be understood that is also applicable.

  As shown in FIG. 1, typical PCR growth curve data is, for example, due to the number of PCR cycles defining the x-axis and an indicator of accumulated polynucleotide growth defining the y-axis. Can be expressed in a two-dimensional coordinate system. Since the use of fluorescent markers is usually the most widely used labeling method, the accumulated growth indicator is the intensity value of the fluorescence, as shown in FIG. . However, it should be understood that other indicators may be used depending on the particular label and / or detection scheme used. Examples of other useful indicators of accumulated signal growth include luminescence intensity, chemiluminescence intensity, bioluminescence intensity, phosphorescence intensity, charge transfer, voltage, current, power, energy, temperature, viscosity, light Scattering, radioactive intensity, reflectance, transmittance, and absorbance are included. The definition of cycle can also include time, process cycle, unit operation cycle, and regeneration cycle.

(Overview of general process)
In accordance with the present invention, one embodiment 100 of a process for determining the transition value of a single sigmoid curve, such as the elbow value or Ct value of a kinetic PCR amplification curve, will be briefly described with reference to FIG. Is possible. In step 110, an experimental data set representing a curve is received or acquired. FIG. 1 shows an example of a plotted PCR data set, where the y-axis and x-axis represent the fluorescence intensity and cycle number of the PCR curve, respectively. In some embodiments, the data set should contain data that is continuous and equidistant along the axis.

  In an exemplary embodiment of the invention, the method includes an input device for entering a data set, such as a keyboard, mouse, and the like; a particular feature of interest in the area of the curve, such as a monitor. A display device for representing points; a processing device such as a CPU necessary for performing each stage of the method; a network interface such as a modem; a data storage device for storing data sets; a computer running on a processor Can be implemented by using a conventional personal computer system, including but not limited to code; and the like. The method can also be performed in a PCR device.

  FIG. 21 shows a system according to the present invention. This figure shows a schematic block diagram showing the relationship between software and hardware resources. The system has a dynamic PCR analysis module that can be placed in a thermocycler and an intelligence module that is part of a computer system. The data set (PCR data set) is transferred from the analysis module to the intelligence module (or vice versa) via a network connection or a direct connection. The data set is processed in accordance with the method shown in FIG. 2 by computer code running on the processor and stored in the intelligence module storage device, and then returned to the analysis module storage device, The changed data can be displayed on the display device. In certain embodiments, an intelligent module may be implemented in the PCR data acquisition device.

  In the case where the process 100 is implemented in an intelligence module (eg, as a processor execution instruction) located in a PCR data acquisition device such as a thermocycler, the data set is collected in real time as the data is collected. Alternatively, it can be stored in a memory unit or buffer and supplied to the intelligence module after the experiment is complete. Similarly, a dataset can be sent to a desktop computer system or other via a network connection (eg, LAN, VPN, intranet, Internet, etc.) or direct connection (eg, USB or other direct wired or wireless connection) to the acquisition device. Can be provided on a separate system such as a computer system or can be provided on a portable medium such as a CD, DVD, floppy disk, or the like. is there. In some aspects, the data set includes data points (or two-dimensional vectors) having coordinate value pairs. In the case of PCR data, the coordinate value pair usually represents the cycle number and the fluorescence intensity value. After the data set is received or acquired in step 110, the end point of the baseline region can be determined by analyzing the data set.

  In step 120, an approximation of the curve is calculated. At this stage, in one embodiment, a double sigmoid function having parameters determined by a Levenberg-Marquardt (LM) regression process or other regression process is used to detect an approximation of the curve representing the data set. Yes. This approximation is referred to as “robust” because outliers or spike points have a minimal impact on the quality of the curve fit. FIG. 13 (which is described below) shows a plot of the received data set and a robust approximation of the data set determined by determining the parameters of the double sigmoid function using the Levenberg-Marquadtr regression process according to the present invention. An example is shown.

  In some aspects, outliers or spike points in the data set are removed or replaced before the data set is processed to determine the endpoint of the baseline region. Spike removal can be performed before or after acquiring the data set at step 110. FIG. 3 shows a process flow for identifying and replacing spike points in a data set representing PCR or other growth curves.

  In step 130, for example, the baseline slope is removed by normalizing the curve using the parameters determined in step 120, as described below. Normalization in this manner allows the Ct value to be determined without having to determine or define the end point of the baseline area of the curve or the stop position of the baseline. As a result, as described below, in step 140, the normalized curve is processed to determine the Ct value.

(LM regression process)
As will be described later, steps 502 to 524 of FIG. 3 also illustrate a process flow (step 120) for approximating the curves of the data set and determining the parameters of the fit function. These parameters can be used for curve normalization (stage 130) (this is the baseline of a data set representing a sigmoid or growth type curve, eg, a PCR curve, according to one embodiment of the invention). Changing or removing the slope). Once the modified data set has been generated with spike points removed or replaced by processing the data set, the fit function is obtained by processing the data set that does not have the modified spike according to steps 502-524. Can be identified.

  In one illustrated embodiment, a robust curve approximation of the data set is calculated by using the Levenberg-Marquardt (LM) method. The LM method is a non-linear regression process, which is an iterative technique that minimizes the distance between the non-linear function and the data set. This is a process that behaves like a combination of a steepest descent process and a Gauss-Newton process, but behaves like a steepest descent process if the current approximation is not well-fitted (which is slow Will behave like a Gauss-Newton process as the current approximation becomes accurate (which is relatively fast, but relative) Convergence with low reliability). The LM regression method is widely used to solve nonlinear regression problems.

  In general, LM regression methods include algorithms that require various inputs and provide outputs. In one aspect, the input includes the data set to be processed, the function used to fit the data, and an initial estimate of the function's parameters or variables. The output includes a set of function parameters that minimize the distance between the function and the data set.

  According to one embodiment, the fit function has the following double sigmoid form:

f (x) = a + bx + c / ((1 + exp- ^{d (xe)} ) (1 + exp- ^{f (xg)} )) (1)

  The reason for choosing this equation as the fit function is based on its flexibility and ability to fit the various curve shapes that a typical PCR curve or other growth curve has. One skilled in the art will appreciate that variations of the aforementioned fit functions and other fit functions can be used as needed.

  The double sigmoid equation (1) has seven parameters, a, b, c, d, e, f, and g. This equation can be decomposed into the sum of a constant, a slope, and a double sigmoid. The double sigmoid itself is a multiplication of two sigmoids. FIG. 4 shows the decomposition of the double sigmoid equation (1). The parameters d, e, f, and g determine the shape of the two sigmoids. To show these effects on the final curve, consider the following single sigmoid.

1 / (1 + exp ^{-d (xe)} ) (2)

  Here, the parameter d determines the “sharpness” of the curve, and the parameter e determines the x value of the inflection point. FIG. 5 shows the influence of the parameter d on the curve and the influence of the parameter e on the position of the inflection point x value. Table 1 below shows the effect of the parameters on the double sigmoid curve.

  In one aspect, the double sigmoid “sharpness” parameters d and f need to be constrained to prevent the curve from acquiring an unrealistic shape. Thus, in one aspect, all iterations at d <−1 or d> 1.1, or f <−1 or f> 1.1 are considered failures. In other aspects, different constraints on the parameters d and f can be used.

  Since the Levenberg-Marquardt algorithm is an iterative algorithm, it usually requires an initial estimate of the function parameters for the fit. The better the initial estimate, the better the approximation and the less likely the algorithm will converge towards a local minimum. Due to the complexity of the double sigmoid function and the various shapes of the PCR curve or other growth curves, only one initial estimate for all parameters can be used to prevent the algorithm often converging to a local minimum. Then it will not be enough. Accordingly, in one aspect, a plurality of (for example, three or more) initial parameter sets are input and the best result is maintained. In one aspect, the majority of parameters are held constant across the plurality of parameter sets in use, and the parameters c, d, and f may differ for each of the plurality of parameter sets. Only. FIG. 6 shows examples of three curve shapes in different parameter sets. The selection of these three parameter sets shows the three possible different shapes of the curve representing the PCR data. It should be understood that more than three parameter sets can be processed to keep the best.

  As shown in FIG. 3, in step 510, initial input parameters for the LM method are identified. These parameters can be input or calculated by the operator. According to one aspect, these parameters are determined or set by steps 502, 504, and 506, as described below.

(Calculation of initial parameter (a))
Parameter (a) is the height of the baseline, and this value is the same across all sets of initial parameters. In one aspect, in step 504, the third smallest y-axis value (eg, fluorescence value) is assigned to parameter (a) from the data set. As a result, stable calculation is possible. Of course, in other embodiments, any other fluorescence value, such as the smallest y-axis value or the second smallest value, can be assigned to parameter (a) as needed.

(Calculation of initial parameter (b))
Parameter (b) is the baseline and plateau slope. This value is the same across all sets of initial parameters. Ideally, no slope should exist, so in one aspect, in step 502, a fixed value of 0.01 is assigned to (b). In other aspects, different values can be assigned to parameter (b), such as values ranging from 0 to about 0.5.

(Calculation of initial parameter (c))
The parameter (c) represents the height of the plateau-the height of the baseline, which is expressed as AFI (Absolute Fluorescence Increase). In one aspect, c = AFI + 2 for the first set of parameters and c = AFI for the last two sets of parameters. This is shown in FIG. 6, where c = AFI for the last two parameter sets and c = AFI + 2 for the first parameter set. This change is due to the shape of the curve modeled by the first set of parameters, which has no plateau.

(Calculation of parameters (d) and (f))
Parameters (d) and (f) define the sharpness of the two sigmoids. In one aspect, in step 502, three typical fixed values are used because there is no way to give an approximation based on the curves in these parameters. It should be understood that other fixed or non-fixed values can be used for parameters (d) and / or (f). These pairs model the common shape of most of the PCR curves encountered. Table 2 below shows the values of (d) and (f) for the different sets of parameters shown in FIG.

(Calculation of parameters (e) and (g))
In step 506, parameters (e) and (g) are determined. Parameters (e) and (g) define two sigmoid inflection points. In one aspect, they all have the same value across all initial parameter sets. Parameters (e) and (g) can have the same or different values. To detect the approximation, in one aspect, the x value of the first point above the average intensity is used (this is, for example, fluorescent but not spiked). A process for determining the values of (e) and (g) according to this aspect is illustrated in FIG. 7 and will be described below.

  Referring to FIG. 7, first, an average of curves (for example, fluorescence intensity) is determined. The first data point above the average is then identified. Then (a) whether the point is not near the start of the curve (eg within the first 5 cycles), (b) the point is near the end of the curve (eg the last 5 And (c) whether the derivative around that point (eg, within the radius of two points around it) does not show a change in sign (showing the change) If so, the point is likely to be a spike and therefore needs to be rejected).

  Table 3 below shows examples of initial parameter values used in FIG. 6 according to one aspect.

  Referring again to FIG. 3, at step 510, after all parameters are set, the LM process 520 is performed using the input data set, function, and parameters. Conventionally, the nonlinear least squares problem is solved by using the Levenberg-Marquardt method. In the conventional LM method, a distance scale defined as the sum of the approximation of a curve and the square of the error between data sets is calculated. However, when minimizing the sum of the squares, the outlier distance is greater than the distance of the non-spike data point, so a large weight is given to the outlier value. Often, inappropriate or undesirable curves will be obtained. Therefore, according to one aspect of the present invention, the distance between the approximation and the data set is calculated by minimizing the sum of absolute errors (as a result, in this case, a large weight is applied to the abnormal value). Is not granted). In this embodiment, the distance between the approximation and the data is given by the following equation.

Distance = Σ | y _data −y _{approximation} | (3)

  As described above, in some aspects, as shown in steps 522 and 524, each of a plurality (eg, three) of sets of initial parameters is entered and processed to retain the best results. This best result is the set of parameters that provides the smallest or smallest distance in equation (3). In one aspect, most of the parameters are held constant across multiple parameter sets, and only c, d, and f can be different for each set of parameters. It should be understood that any number of initial parameter sets can be used.

  FIG. 8 shows the process flow of the LM process 520 for a set of parameters according to the present invention. As described above, the Levenberg-Marquardt method can operate like a steepest descent process or a Gauss-Newton process. This behavior depends on the attenuation coefficient λ. The larger λ, the more the Levenberg-Marquardt algorithm will behave like a steepest descent process. On the other hand, as λ is smaller, the Levenberg-Marquardt algorithm behaves like a Gauss-Newton process. In one aspect, λ starts at 0.001. It should be understood that λ can start from any other value, such as from about 0.000001 to about 1.0.

  As described above, the Levenberg-Marquardt method is an iterative technique. According to one aspect, the following steps are performed in each iteration, as shown in FIG.

  1. The preceding approximate Hessian matrix (H) is calculated.

2. The preceding approximate transposed Jacobian matrix (J ^T ) is calculated.

  3. The preceding approximate distance vector (d) is calculated.

  4). With the current attenuation coefficient λ, the diagonal term of the Hessian is increased as follows.

H _aug = Hλ (4)

  5. Solve the augmented equation as follows.

H _aug x = J ^T d (5)

  6). Add the increased equation solution x to the function parameters.

  7). Calculate new approximations and distances between curves.

  8). If the distance with the new parameter set is less than the distance with the previous parameter set, then the iteration is considered successful, and keeps or saves the new set of parameters. Reduce λ (eg, by a factor of 10). If the distance with the new parameter set is greater than the distance with the previous parameter set, then the iteration is considered a failure, discards the new set of parameters, and sets the damping factor λ to Increase (for example, by a factor of 10).

  In one aspect, the LM process of FIG. 8 is repeated until one of the following criteria is achieved:

  1. Already executed for a predetermined number of iterations N. This first criterion prevents the algorithm from being repeated indefinitely. For example, in one aspect shown in FIG. 10, the default iteration value N is 100. If the algorithm can converge, 100 iterations should be sufficient for the algorithm to converge. In general, N can range from less than 10 to more than 100.

  2. The difference in distance between two successful iterations is below a threshold (eg, 0.0001). If the difference becomes very small, it is pointless to continue the iteration because the desired accuracy has already been achieved and the solution will not be much better.

3. The attenuation coefficient λ exceeds a predetermined value (for example, exceeds 10 ²⁰ ). When λ becomes very large, the algorithm does not converge better than the current solution, and therefore it does not make sense to continue iterating. In general, this specified value may be significantly below or above 10 ²⁰ .

(Normalization)
After the parameter determination is complete, in one embodiment, the curve is normalized using one or more of the determined parameters (step 130). For example, in one aspect, a curve can be normalized or adjusted to have a zero baseline slope by subtracting the linear growth portion of the curve. Mathematically, this is expressed as:

dataNew (BLS) = data− (a + bx) (6)

  Here, dataNew (BLS) is a normalized signal after baseline subtraction, which is, for example, a data set (data) from which linear growth or baseline slope has been subtracted or removed. The values of parameters a and b are values determined by curve regression using the LM equation, and x is the number of cycles. Therefore, for all data values along the x-axis, a data curve having a zero baseline slope is generated by subtracting the constant a and slope b multiplied by the x value from the data. In certain aspects, spike points are removed from the data set before applying the LM regression process to the data set to determine the normalization parameters.

  In another aspect, the curve can be normalized or adjusted to have a zero slope according to the following equation:

dataNew (BLSD) = (data- (a + bx)) / a (7a)

  Where dataNew (BLSD) is the normalized signal after baseline subtraction with division, which subtracts or removes the linear growth or baseline slope, for example, and divides the result by a Data set (data). The values of parameters a and b are values determined by curve regression using the LM equation, and x is the number of cycles. Thus, for all values along the x axis, the constant a and the slope b multiplied by the x value are subtracted from the data, and the result is divided by the value of the parameter a, thereby having a baseline slope of zero. A data curve is generated. In one aspect, the equation (7a) is effective when the parameter “a” ≧ 1, and when “a” <1, the following equation is used.

dataNew (BLSD) = data- (a + bx) (7b)

  In certain aspects, spike points are removed from the data set before applying the LM regression process to the data set to determine the normalization parameters.

  In yet another aspect, the curve can be normalized or adjusted according to the following equation:

dataNew (BLD) = data / a (8a)

  Here, dataNew (BLD) is a normalized signal after baseline division, which is, for example, a data set (data) divided by parameter a. The values of parameters a and b are values determined by curve regression using the LM equation, and x is the number of cycles. In one aspect, equation (8a) is valid when parameter “a” ≧ 1, and when “a” <1, the following equation is used.

dataNew (BLD) = data + (1-a) (8b)

  In certain aspects, spike points are removed from the data set before applying the LM regression process to the data set to determine the normalization parameters.

  In yet another aspect, the curve can be normalized or adjusted according to the following equation:

dataNew (PGT) = (data- (a + bx)) / c (9a)

  In this case, dataNew (PGT) is the normalized signal after baseline subtraction with division, which subtracts or removes the linear growth or baseline slope, for example, and divides the result by c Data set (data). The values for a, b, and c are values determined by curve regression using the LM equation, and x is the number of cycles. Thus, for all data values along the x axis, the constant a and the slope b multiplied by the x value are divided from the data, and the result is divided by the value of the parameter c, so that a zero baseline slope is obtained. A data curve is generated. In one aspect, the equation (9a) is effective when the parameter “c” ≧ 1, and the following equation is used when the parameters “c” <1 and “c” ≧ 0.

dataNew (PGT) = data- (a + bx) (9b)

  In one aspect, spike points are removed from the data set before applying the LM regression process to the data set to determine the normalization parameters.

  Those skilled in the art will appreciate that other normalization formulas can be used to normalize and / or modify the baseline by using parameters determined by Levenberg-Marquardt or other regression processes. Will.

(Decision of curvature)
The Ct value can be determined after normalization of the curve using equations (6), (7), (8), (9), or any other normalization equation is complete. In one aspect, the curvature determination process or method is transformed into a normalized curve, as described with reference to FIG. 9, which shows a process flow for determining the elbow or Ct value of a kinetic PCR curve. Applicable. In step 910, a data set is acquired. In cases where the decision process is being performed within an intelligence module (eg, as a processor execution instruction) located within a PCR data acquisition device such as a thermocycler, the data set is transferred to the intelligence module as data is collected. It can be supplied in real time, or it can be stored in a memory unit or buffer and supplied to the module after the experiment is complete. Similarly, the data set can be sent to the acquisition device via a network connection (eg, LAN, VPN, intranet, internet, etc.) or a direct connection (eg, USB or other direct wired or wireless connection), etc. It can be provided on a separate system, or it can be provided on a portable medium such as a CD, DVD, floppy disk, or the like.

  After the receipt or acquisition of the data set is complete, at step 920, an approximation to the curve is determined. At this stage, in one embodiment, an approximation of the curve representing the data set is detected by using a double sigmoid function having parameters determined by the Levenberg-Marquardt regression process. It is also possible to remove spike points from the data set prior to step 920, as described with reference to FIG. For example, the data set acquired in step 910 may be a data set from which spikes have already been removed. In step 930, the curve is normalized. In some embodiments, the curve is normalized using any of Equations (6), (7), (8), or (9) described above. For example, the baseline can be set to zero slope by subtracting the baseline slope according to equation (6) above using the double sigmoid parameter determined in step 920. In step 940, the process is applied to the normalized curve to determine the curvature at points along the normalized curve. A plot of curvature versus cycle number can be returned and / or displayed. The point of maximum curvature corresponds to the elbow value or Ct value. In step 950, the results are returned to, for example, the system that performed the analysis or a separate system that requested the analysis. In step 960, the Ct value is displayed. It is also possible to display additional data such as the entire data set or curve approximation. The graphical display can be rendered by a display device, such as a monitor screen or printer, coupled to the system that performed the analysis of FIG. 9, or data is provided to a separate system for rendering on the display device. It is also possible to do.

  According to one embodiment, the maximum curvature is determined to obtain the Ct value of this curve. In one aspect, the curvature at some or all points on the normalized curve is determined. A plot of curvature versus number of cycles can be displayed. The curvature of the curve is given by:

  Consider a circle of radius a given by:

  The curvature of equation (11) is kappa (x) = − (1 / a). Therefore, the radius of curvature is equal to the negative reciprocal of the curvature. Since the radius of the circle is constant, its curvature is given by-(1 / a). Now consider FIG. 10b. This is a plot of the curvature of the fit of the PCR data set of FIG. 10a. The Ct value is considered to occur at the position of maximum curvature, which occurs at cycle number Ct = 21.84. This Ct value corresponds well with the PCR growth curve shown in FIG. 10a.

  The radius of curvature at the maximum curvature (which corresponds to a Ct value of 21.84) is Radius = 1 / 0.2818 = 3.55 cycles. FIG. 11 shows the circle of this radius superimposed on the PCR growth curve of FIG. 10a. As shown in FIG. 11, the circle with the radius corresponding to the maximum curvature represents the maximum circle that can be superimposed on the starting point of the growth region of the curve while being tangent to the curve. A curve with a small (maximum) radius of curvature will have a steep growth curve, and a curve with a large (maximum) radius of curvature will have a gentle growth curve. If the radius of curvature is extremely large, this indicates a curve with no obvious signal.

  The first and second derivatives of the double sigmoid of the formula (1) necessary for calculating the curvature are shown below.

(First derivative)

(Second derivative)

  FIG. 12a shows an example of raw data of the growth curve. By applying the double sigmoid / LM method to the raw data plot shown in FIG. 12b, the values of the seven parameters of equation (1) shown in Table 4 below are obtained.

  FIG. 13 shows a double sigmoid fit to the data shown in FIG. 12, which shows a very accurate assessment of the data points. Then, by normalizing these data according to equation (6) (baseline subtraction), the graph shown in FIG. 14 is obtained. The solid line shown in FIG. 14 is the result of applying the double sigmoid / LM of equation (1) to this data set normalized according to equation (6). FIG. 15 shows a plot of curvature versus number of cycles in the normalized curve of FIG. A curvature of 0.1378, which is the maximum value of curvature, occurs at a cycle number of 34.42. Therefore, based on the number of cycles at this maximum curvature, Ct = 34.42 and the radius of curvature = 1 / 0.1378 = 7.25. FIG. 16 shows a superposition of a circle having this radius of curvature and a normalized data set.

  An example of a “slow-grower” data set is shown in FIG. Normalization using double sigmoid fit and baseline subtraction (Equation (6)) on this data set yields the fit results shown in FIG. The corresponding curvature plot is shown in FIG. Curvature = 0.00109274, which is the maximum curvature, occurs at a cycle number of 25.90, which corresponds to a radius of curvature = 915. This large radius of curvature indicates that this is a slow growth data set.

  As another example, consider the set of PCR growth curves shown in FIG. A comparison of Ct values obtained using the existing method (“threshold”) versus baseline subtraction with division (BLSD—equation (7)) is shown in Table 5 below. It is shown.

  Table 5 shows that according to the curvature method for calculating the Ct value (in this case, after normalization by BLSD), a smaller Cv (Coefficient of Variation) is obtained compared to the existing threshold method. It is shown that. The radius of curvature (ROC) calculated by the curvature method provides a simple method for discriminating between a linear curve and a true growth curve.

(Conclusion)
According to one aspect of the invention, a computer-implemented method is provided for determining a point at the end of a baseline of a growth curve. The method typically includes receiving a data set representing a growth curve, the data set including a plurality of data points each having a pair of coordinate values, and a Levenberg-Marquardt (LM). ) Determining the parameters of the function by applying a regression process to the double sigmoid function to calculate an approximation of the curve that fits the data set. The method typically involves normalizing the curve using the determined parameters to generate a normalized curve, and processing the normalized curve to determine the coordinate value of the point of maximum curvature. And this point represents the end point of the baseline region of the growth curve. In one aspect, the normalized curve is processed to determine the curvature at some or all points along the curve. In certain embodiments, the method further includes displaying a plot of the curvature of the normalized curve.

  In one aspect, the data set represents a growth curve of a kinetic polymerase chain reaction (PCR) process, a bacterial process, an enzymatic process, or a binding process. In certain embodiments, the accumulation of amplified polynucleotide comprises fluorescence intensity values, luminescence intensity values, chemiluminescence intensity values, phosphorescence intensity values, charge transfer values, bioluminescence intensity values, or absorption values. Represented by one of them.

  In another example, the curve represents a growth curve for a kinetic polymerase chain reaction (PCR) process, and the point at the end of the baseline region is the elbow or cycle threshold ( Ct), and in this case, the point having the maximum curvature represents the Ct value. In some embodiments, the pair of coordinate values represents the accumulation of amplified polynucleotide and the number of cycles. In another aspect, the method further includes returning a Ct value. In yet another aspect, the method further includes displaying a Ct value.

  In some embodiments, normalizing includes subtracting a linearly growing portion of the curve. In some aspects, the received data set includes a data set that has been processed to remove one or more outliers or spike points.

In one embodiment, the double sigmoid function has the form a + bx + c / ((1 + exp ^{−d (xe)} ) (1 + exp ^{−f (xg)} )), and the calculation step includes parameters a, b, including iteratively determining one or more of c, d, e, f, and g. In a particular embodiment, at least parameters a and b are determined, and the normalization step includes subtracting the linear growth portion a + bx from the curve. In another specific embodiment, at least the parameter a is determined, and the normalizing step includes dividing the curve by the parameter a. In yet another specific embodiment, at least parameters a and b are determined and normalizing includes subtracting the linear growth portion a + bx from the curve and dividing the result by parameter a.

  In another embodiment, at least parameters a, b, and c are determined, and the normalization step includes subtracting the linear growth portion a + bx from the curve and dividing the result by the parameter c. Yes.

  In accordance with another aspect of the invention, a computer readable medium is provided that includes code for controlling a processor to determine a point at the end of a baseline region of a growth curve. This code typically receives a data set that represents a growth curve (this data set includes multiple data points, each with a pair of coordinate values), and doubles the Levenberg-Marquardt (LM) regression process. It includes instructions for determining the parameters of this function by calculating an approximation of the curve that is applied to the function to fit the data set. The code also typically uses the determined parameters to normalize the curve to generate a normalized curve, and process the normalized curve to determine the coordinate value of the point of maximum curvature. This point also represents the end point of the baseline region of the growth curve.

  In certain embodiments, the computer readable medium further includes instructions for displaying a normalized curve curvature plot. In certain aspects, the code further includes instructions for returning or displaying the coordinate value of the point at the end of the baseline region.

  In one aspect, the data set represents a growth curve of a kinetic polymerase chain reaction (PCR) process, a bacterial process, an enzymatic process, or a binding process. In another aspect, the computer readable medium further includes instructions for displaying the Ct value. In yet another embodiment, the computer readable medium further includes instructions for returning a Ct value.

  In certain embodiments, the data set represents a growth curve of a kinetic polymerase chain reaction (PCR) process, and the point at the end of the baseline region represents the elbow or cycle threshold (Ct) of the growth curve. ing. In another embodiment, the pair of coordinate values represents the accumulation of amplified polynucleotide and the number of cycles. In yet another embodiment, the accumulation of amplified polynucleotide comprises fluorescence intensity values, luminescence intensity values, chemiluminescence intensity values, phosphorescence intensity values, charge transfer values, bioluminescence intensity values, or absorption values. Represented by one of the

  In one aspect, the instructions for processing include instructions for determining curvature at some or all points along the normalized curve. In certain aspects, the instructions for normalization include instructions for subtracting the linear growth portion from the data set. In another aspect, the curve is an amplification curve of a kinetic polymerase chain reaction (PCR) process, and the point at the end of the baseline region represents the elbow or cycle threshold (Ct) of the kinetic PCR curve. The point having the maximum curvature represents the Ct value.

In a particular embodiment, the double sigmoid function has the form a + bx + c / ((1 + exp ^{−d (xe)} ) (1 + exp ^{−f (xg)} )), and the instruction to calculate is Instructions are included for iteratively determining one or more of the parameters a, b, c, d, e, f, and g.

  In some aspects, at least parameters a and b are determined, and the instructions for normalization include instructions for subtracting the linear growth portion a + bx from the curve. In another aspect, at least parameters a, b, and c are determined, and the instructions for normalization include instructions for subtracting the linear growth portion a + bx from the curve and dividing the result by parameter c. Yes. In yet another aspect, at least the parameter a is determined, and the instruction for normalizing includes an instruction for dividing the curve by the parameter a. In yet another aspect, at least parameters a and b are determined, and the instructions for normalization include instructions for subtracting the linear growth portion a + bx from the curve and dividing the result by parameter a.

  According to yet another aspect of the invention, a kinetic polymerase chain reaction (PCR) system is provided. The system typically includes a kinetic PCR analysis module that generates a PCR data set that represents a kinetic PCR amplification curve (this data set includes a plurality of data points, each having a pair of coordinate values, And the data set includes a data point within a region of interest that includes a cycle threshold (Ct)) and an intelligence module adapted to process the PCR data set to determine a Ct value. It is out. The intelligence module typically determines the parameters of the function by applying a Levenberg-Marquardt (LM) regression process to the double sigmoid function to calculate a curve fit that fits the data set, and uses the determined parameters. Normalizing the curve to generate a normalized curve and processing the normalized curve to determine the coordinate value of the point of maximum curvature (this point is the growth curve cycle threshold (Ct)) Represent the PCR data set.

In one embodiment, the double sigmoid function has the form a + bx + c / ((1 + exp ^{−d (xe)} ) (1 + exp ^{−f (xg)} )), and the calculation step includes parameters a, b, including iteratively determining one or more of c, d, e, f, and g. In a particular embodiment, at least parameters a and b are determined, and the normalization step includes subtracting the linear growth portion a + bx from the curve. In another specific embodiment, the parameters a and b are determined and normalizing includes subtracting the linear growth portion a + bx from the curve and dividing the result by the parameter a. In yet another specific embodiment, at least the parameter a is determined and normalizing includes dividing the curve by the parameter a.

  In another certain aspect, at least parameters a, b, and c are determined and normalizing comprises subtracting the linear growth portion a + bx from the curve and dividing the result by parameter c. It is out.

  In one aspect, the normalized curve is processed to determine the curvature at some or all points along the curve, where the point with the maximum curvature represents the Ct value. . In certain embodiments, the intelligence module is further adapted to render a display of a normalized curve curvature plot.

  In some embodiments, normalizing includes subtracting the linear growth portion from the data set. In certain embodiments, the curve fit is normalized by subtracting the linear growth portion a + bx from the curve.

  In certain embodiments, the intelligence module is further adapted to return Ct values. In another embodiment, the intelligence module is further adapted to display Ct values.

  In another embodiment, the coordinate value pair represents the accumulation of amplified polynucleotide and the cycle number. In certain embodiments, the accumulation of amplified polynucleotide comprises fluorescence intensity values, luminescence intensity values, chemiluminescence intensity values, phosphorescence intensity values, charge transfer values, chemiluminescence intensity values, or absorption values. Represented by one of them.

  In certain embodiments, the kinetic PCR analysis module resides in a kinetic thermocycler device, and the intelligence module includes a processor communicatively coupled to the analysis module. In another specific embodiment, the intelligence module includes a processor residing in a computer system coupled to the analysis module by either a network connection or a direct connection.

  It should be understood that the Ct determination process, including the curve fitting and curvature determination process, can be implemented as computer code running on a processor of a computer system. This code includes instructions for controlling the processor to implement various aspects and steps of the Ct determination process. This code is usually stored on a hard disk, RAM, or a portable medium such as a CD or DVD. Similarly, these processes may be performed in a PCR device such as a thermocycler that includes processor execution instructions stored in a memory unit coupled to the processor. The code containing such instructions can be downloaded to the PCR device memory via a network connection or a direct connection to the code source, or by using a portable medium, as is well known. .

  For those skilled in the art, the elbow determination process of the present invention can be implemented in various programming languages such as C, C ++, C #, Fortran, VisualBasic, and prepackaged routines, functions useful for data visualization and analysis. It will be understood that it can be coded using an application such as Mathematica, which provides the procedure. A further example of the latter is MATLAB®.

  Although the invention has been described above by way of example in terms of specific embodiments, it should be understood that the invention is not limited to the disclosed examples. On the contrary, the present invention is intended to cover various modifications and similar arrangements that will be apparent to those skilled in the art. Accordingly, the scope of the appended claims should be accorded the broadest interpretation so as to encompass all such modifications and similar arrangements.

Figure 2 shows an example of a typical PCR growth curve plotted as fluorescence intensity versus cycle number. FIG. 5 shows a process flow for determining the end point of the baseline region of a growth curve (or Ct value of a PCR curve). Fig. 4 shows a detailed process flow of a spike identification and replacement process according to an embodiment of the invention. Fig. 4 shows a detailed process flow of a spike identification and replacement process according to an embodiment of the invention. Fig. 4 shows a detailed process flow of a spike identification and replacement process according to an embodiment of the invention. Fig. 4 shows a detailed process flow of a spike identification and replacement process according to an embodiment of the invention. Fig. 4 shows a detailed process flow of a spike identification and replacement process according to an embodiment of the invention. Fig. 4 shows a detailed process flow of a spike identification and replacement process according to an embodiment of the invention. Fig. 4 illustrates a double sigmoid decomposition including parameters a to g. The influence of parameter (d) on the curve and the influence of (e) on the position of the inflection point x value are shown. All the curves in FIG. 5 have the same parameters except for parameter (d). An example of three curve shapes in different parameter sets is shown. FIG. 6 illustrates a process for determining values of double sigmoid parameters (e) and (g) according to one aspect. FIG. 6 illustrates a process for determining values of double sigmoid parameters (e) and (g) according to one aspect. Fig. 4 shows the process flow of a Levenberg-Marquardt regression process with an initial set of parameters. Fig. 4 shows the process flow of a Levenberg-Marquardt regression process with an initial set of parameters. Fig. 4 shows the process flow of a Levenberg-Marquardt regression process with an initial set of parameters. Fig. 4 shows the process flow of a Levenberg-Marquardt regression process with an initial set of parameters. FIG. 4 illustrates a more detailed process flow for determining the elbow value of a PCR process according to an embodiment. FIG. 10a shows a typical growth curve fitted to experimental data using double sigmoid. FIG. 10b shows a plot of the curvature of the double sigmoid curve of FIG. 10a. 10a shows a circle tangent to the point of maximum curvature superimposed on the growth curve of FIG. 10a. An example of a growth curve data set is shown. Fig. 12a shows a plot of the data set of Fig. 12a. FIG. 13 shows a double sigmoid fit for the data set of FIG. FIG. 13 shows the data set of FIG. 12 (FIG. 13) (and the double sigmoid fit) after normalization using the baseline subtraction method of equation (6). FIG. 15 shows a plot of curvature versus cycle number for the normalized data set of FIG. FIG. 15 shows the superposition of the circle with the maximum radius of curvature and the normalized data set of FIG. An example of a “slow-grower” data set is shown. Figure 18 shows the data set and double sigmoid fit of Figure 17 after normalization using the baseline subtraction method of equation (6). FIG. 19 shows a plot of curvature versus number of cycles in the normalized data set of FIG. Figure 5 shows a plot of a set of PCR growth curves including repeated runs and negative samples. 1 shows a schematic block diagram illustrating the relationship between software and hardware resources.

Claims (22)

In a computer implemented method for determining a point at the end of a baseline region of a growth curve,
Receiving a data set representing a growth curve, the data set comprising a plurality of data points each having a pair of coordinate values;
Determining a parameter of the function by applying a Levenberg-Marquardt (LM) regression process to a double sigmoid function to calculate an approximation of a curve that fits the data set;
Normalizing the curve using the determined parameters to generate a normalized curve;
Determining a point of maximum curvature representing the end point of the baseline region of the growth curve by processing the normalized curve;
Having a method.   The method of claim 1, wherein normalizing includes subtracting a linear growth portion from the data set.   The method of claim 1, wherein processing comprises determining the curvature at some or all points along the normalized curve.   The method of claim 3, further comprising displaying a plot of the curvature of the normalized curve. The double sigmoid function has the form a + bx + c / ((1 + exp ^{−d (xe)} ) (1 + exp ^{−f (xg)} )), and the step of calculating includes the parameters a, b, c of the function 2. The method of claim 1, comprising iteratively determining one or more of d, e, f, and g.   6. The method of claim 5, wherein at least the parameters a and b are determined and normalizing comprises subtracting a linear growth portion a + bx from the curve.   6. The method of claim 5, wherein at least the parameter a is determined and normalizing comprises dividing the curve by the parameter a.   6. The method of claim 5, wherein at least the parameters a and b are determined and normalizing comprises subtracting the linear growth portion a + bx from the curve and dividing the result by the parameter a. Method.   6. The step of determining and normalizing at least the parameters a, b, and c includes subtracting the linear growth portion a + bx from the curve and dividing the result by the parameter c. the method of.   The data set represents a kinetic polymerase chain reaction (PCR) process growth curve, and the point at the end of the baseline region represents the elbow or cycle threshold (Ct) of the growth curve. The method of claim 1.   The method of claim 10, further comprising returning the Ct value.   The method of claim 10, further comprising displaying the Ct value. A computer readable medium including code for controlling a processor to determine a point at the end of a baseline region of a growth curve,
The code is
Receiving a data set representing a growth curve, the data set comprising a plurality of data points each having a pair of coordinate values;
Determining a parameter of the function by applying a Levenberg-Marquardt (LM) regression process to a double sigmoid function to calculate an approximation of a curve that fits the data set;
Normalizing the curve using the determined parameters to generate a normalized curve;
Determining a point of maximum curvature representing the end point of the baseline region of the growth curve by processing the normalized curve;
A computer readable medium containing instructions for executing   The computer-readable medium of claim 13, wherein the instructions for processing include instructions for determining the curvature at some or all points along the normalized curve. The double sigmoid function has a form of a + bx + c / ((1 + exp ^{−d (xe)} ) (1 + exp ^{−f (xg)} )), and the instruction to calculate the parameter a of the function The computer-readable medium of claim 13, comprising instructions for iteratively determining one or more of b, c, d, e, f, and g.   The computer-readable medium of claim 13, wherein the code further comprises instructions for returning or displaying the coordinate value of the point at the end point of the baseline region. In a kinetic polymerase chain reaction (PCR) system,
A PCR data set representing a kinetic PCR amplification curve, each including a plurality of data points having coordinate value pairs, and including data points within a region of interest including a cycle threshold (Ct) A kinetic PCR analysis module for generating a running data set;
Determining a parameter of the function by applying a Levenberg-Marquardt (LM) regression process to a double sigmoid function to calculate an approximation of a curve that fits the data set;
Normalizing the curve using the determined parameters to generate a normalized curve;
Determining a point of maximum curvature representing the cycle threshold (Ct) of the growth curve by processing the normalized curve;
An intelligence module adapted to process the PCR data set to determine the Ct value;
Having a system. The double sigmoid function has the form a + bx + c / ((1 + exp ^{−d (xe)} ) (1 + exp ^{−f (xg)} )), and the step of calculating includes the parameters a, b, c of the function 18. The system of claim 17, comprising iteratively determining one or more of d, e, f, and g.   The system of claim 17, wherein processing comprises determining the curvature at some or all points along the normalized curve.   The system of claim 19, wherein the intelligence module is adapted to render a display of the curvature plot of the normalized curve.   The system of claim 17, wherein the kinetic PCR analysis module is present in a kinetic thermocycler device, and the intelligence module includes a processor communicatively coupled to the analysis module.   The system of claim 17, wherein the intelligence module includes a processor residing in a computer system coupled to the analysis module by either a network connection or a direct connection.

Images