An automatic calculation method for elbow point CT value in quantitative polymerase chain reaction

Abstract

The present invention discloses an automatic calculation method for the elbow point CT value of quantitative polymerase chain reaction, which relates to the field of bioinformatics technology and mainly includes the following steps: performing amplification curve typing under linear determination according to the original fluorescence value, screening approximate elbow point positions based on the curve change characteristics of the exponential amplification period of the amplification curve, and screening approximate knee point positions based on the curve change characteristics of the saturation plateau period of the amplification curve; fitting the original fluorescence value amplification curve between the approximate elbow point and the approximate knee point, and based on the optimally fitted equation as the linear region equation; fitting the original fluorescence value amplification curve before the approximate elbow point, and based on the optimally fitted equation as the baseline region equation; solving the intersection point based on the linear region equation and the baseline region equation, and taking this intersection point as the target elbow point CT value. The result obtained by the present invention fully reflects the characteristics of the original amplification curve, and the obtained quantitative result is more stable.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and particularly to an automatic calculation method for the elbow point CT value of quantitative polymerase chain reaction. Background Art

[0002] Real-time fluorescence quantitative polymerase chain reaction (PCR) technology has become one of the main means of molecular research due to its advantages of being real-time, sensitive, and specific, and is widely used in molecular biology, medical research, and clinical practice. However, ordinary real-time quantitative PCR is limited by the fluorescence detection channel, and only a single target sequence can be analyzed in a single reaction, which limits its use in scientific research. The multiplex quantitative PCR method based on multi-color fluorescence detection and analysis technology can detect multiple target sequences in one reaction by labeling different fluorescent dyes in one reaction system.

[0003] Currently, most of the instruments for quantitative calculation adopt the method of manual threshold, which requires manual judgment and removal of various abnormal conditions. The traditional method for automatically obtaining the CT value is to take 8-10 times the variance of a section of baseline values before amplification as the reference baseline for finding the elbow position (elbow point). However, the problem with this algorithm is obvious that it will be affected by the instability of the previous baseline. Especially in the first 5 cycles, the interaction between the initial temperature of the reagent, the internal temperature of the instrument, and the characteristics of the reagent itself leads to large fluctuations during this period, which often extends to the following cycles, and often shows a downward or upward slope.

[0004] Generally speaking, the methods for CT value calculation include mathematical algorithms such as expectation maximization, nearest neighbor analysis, basic model parameterization, Bayesian estimation, or combinations thereof. These mathematical algorithms are the methods actually used in many devices, but in fact, no method can surpass the manual threshold method. The reason is that these methods start more from mathematics, either masking the original state of the signal or sacrificing many signal characteristics during the amplification process. Therefore, although the qPCR technology has been applied for many years, it is still a difficult point for the real-time fluorescence quantitative PCR technology to automatically obtain accurate quantification. Summary of the Invention

[0005] In order to achieve the automatic measurement of real-time fluorescence quantitative PCR, the present invention proposes an automatic calculation method for the elbow point CT value of quantitative polymerase chain reaction, including the steps of:

[0006] S1: Obtain the original fluorescence value amplification curve at a preset number of cycles during the quantitative polymerase chain reaction process;

[0007] S2: Make a preliminary judgment on the amplification curve according to the fluorescence values of the previous preset number of cycles, and enter the next step after judging that there is an effective amplification curve, otherwise end;

[0008] S3: Perform amplification curve typing under linear determination based on the original fluorescence values, screen for approximate elbow point positions based on the curve change characteristics in the exponential amplification phase of the amplification curve, and screen for approximate knee point positions based on the curve change characteristics in the saturation plateau phase of the amplification curve;

[0009] S4: Fit the original fluorescence value amplification curve between the approximate elbow point and the approximate knee point, and the equation after the optimal fit is the linear region equation;

[0010] S5: Fit the original fluorescence value amplification curve before the approximate elbow point, and the equation after the optimal fit is the baseline region equation;

[0011] S6: Solve for the intersection point based on the linear region equation and the baseline region equation, and take this intersection point as the target elbow point CT value.

[0012] Further, in the step S2, the fluorescence values of the preset number of cycles before fitting are made into a straight line, and the difference is compared with the corresponding original fluorescence amplification curve. When the difference is continuously greater than zero starting from any point, it is determined that there is an effective amplification curve.

[0013] Further, in the step S3, the amplification curve typing includes a baseline region, an exponential amplification region, a linear growth region, a logarithmic growth region, and a saturation plateau region.

[0014] Further, in the step S3, the start and end points with a stable fluorescence value change trend between cycle positions are searched for to perform the typing of the amplification curve. Among them, the fluorescence value change trend between cycle positions is determined by the following formula:

[0015]

[0016] In the formula, is the fluorescence value change trend between cycle positions, c is a ladder constant greater than or equal to 1, is the fluorescence value at the nth cycle, and the start and end points are obtained by searching for the corresponding cycle positions of .

[0017] Further, in the step S3, the curve change characteristics in the exponential amplification phase of the amplification curve are that it enters the exponential growth region with exponential growth of the fluorescence value from the low fluorescence value state in the baseline region. Among them, the fluorescence value growth speed trend of the elbow point is determined by the following formula:

[0018]

[0019] In the formula, is the fluorescence value growth speed trend of the elbow point, and by searching in the original fluorescence value amplification curve for the starting low fluorescence value and then within the range of exponential growth Screen the approximate elbow point positions corresponding to the cyclic positions.

[0020] Further, in the step S3, the curve change characteristics of the saturation plateau of the amplification curve are that the fluorescence value in the linear region linearly increases, then turns into the logarithmic growth region and finally enters the saturation plateau region with stagnant growth. Among them, the growth rate trend of the fluorescence value at the knee point is determined by the following formula:

[0021]

[0022] In the formula, M is the maximum fluorescence value, N is the system noise baseline, is the fluorescence value at the nth cycle after correction, is the growth rate trend of the fluorescence value at the knee point, and the approximate knee point positions are screened by searching for the cyclic positions corresponding to the logarithmic growth of the fluorescence value and finally entering the range of stagnant growth in the original fluorescence value amplification curve Screen the approximate elbow point positions corresponding to the cyclic positions.

[0023] Further, in the step S4, the equation of the fitted linear region is expressed as:

[0024]

[0025] In the step S5, the equation of the fitted baseline region is expressed as:

[0026]

[0027] In the formula, x is the number of cycles, Y is the fluorescence value, is the fitting parameter under the optimal fitting of the linear region equation, is the fitting parameter under the optimal fitting of the baseline region equation, and B is the user-defined threshold.

[0028] Compared with the prior art, the present invention has at least the following beneficial effects:

[0029] The automatic elbow point CT value calculation method for quantitative polymerase chain reaction according to the present invention uses the characteristics of line segments for segmentation, then performs fitting after segmentation, and uses the fitting line characteristics to calculate CT with a simple and effective model. Calculating the CT value using the fitting line characteristics of the segmented segments, its result fully reflects the characteristics of the original amplification curve, and it has better compliance characteristics and more realistic effects for low-concentration samples such as weak positives. At the same time, the CV of the relative quantitative result is more stable. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 is a step diagram of an automatic elbow point CT value calculation method for quantitative polymerase chain reaction;

[0031] Figure 2 is a schematic diagram of PCR amplification curve and automatic calculation principle analysis. Detailed implementation manners

[0032] The following are specific embodiments of the present invention. In combination with the accompanying drawings, the technical solutions of the present invention will be further described, but the present invention is not limited to these embodiments.

[0033] During the PCR amplification process, signals such as inclination and protrusion appear because although the instrument itself has a thermal cover to avoid the influence of condensation, liquid evaporation and condensation still occur during rapid temperature change. Therefore, it can be seen that the amplified liquid condenses on an annular tube wall adjacent to the bottom liquid. Whether fluorescence is collected from above or below, due to the liquid state and position, it will cause the continuous sinking of the fluorescence curve during the amplification process. In response to this phenomenon, ABI7500 uses a ROX correction strategy to correct it. However, due to the change in the liquid volume, the heat capacity decreases, resulting in the actual temperature in the later stage of the amplification segment being slightly higher than that in the earlier stage, thereby causing a change in the enzyme efficiency and making the amplified fluorescence signal show an upward or downward change, resulting in a very complex state of the actual amplification curve.

[0034] Currently, a relatively good method for automatic CT value calculation is to fit the original data using a four-parameter or five-parameter curve fitting method, or directly use the maximum second derivative to find the position of the curve change, and it also includes a method of first fitting to form new data and then using the maximum second derivative for calculation. Whether the fitted data can reflect the characteristics of the original amplification curve, the original amplification signal is discarded, and the signal generated by the fitting formula is used instead, which makes the user unable to see what the actual physical signal characteristics are, resulting in the user not knowing the actual effect of their reagent during the development of the kit, and instead there are potential problems.

[0035] If you want to use this method without distortion, the fitting line must be able to remove adverse conditions such as inclination, protrusion, and starting mutations while trying to restore and display the basic characteristics of the original signal as much as possible. However, for the four-parameter and five-parameter S curve fitting, the error value with the original signal is very large. Although this method can obtain more stable CT values for large signals, once the sensitivity is measured, there will be many abnormal data. For the problems existing in the prior art, as Figure 1 shown, the present invention proposes an automatic calculation method for the elbow point CT value of quantitative polymerase chain reaction, including the steps:

[0036] S1: Obtain the original fluorescence value amplification curve at a preset number of cycles during the quantitative polymerase chain reaction process;

[0037] S2: Make a preliminary judgment on the amplification curve according to the fluorescence values of the previous preset number of cycles, and enter the next step after judging that there is an effective amplification curve, otherwise end;

[0038] S3: Perform amplification curve typing under linear determination based on the original fluorescence values, screen for approximate elbow points based on the curve change characteristics during the exponential amplification period of the amplification curve, and screen for approximate knee points based on the curve change characteristics during the saturation plateau period of the amplification curve;

[0039] S4: Fit the original fluorescence value amplification curve between the approximate elbow point and the approximate knee point, and the equation after the optimal fit is the linear region equation;

[0040] S6: Fit the original fluorescence value amplification curve before the approximate elbow point, and the equation after the optimal fit is the baseline region equation;

[0041] S6: Solve for the intersection point based on the linear region equation and the baseline region equation, and use this intersection point as the target elbow point CT value.

[0042] As Figure 2 shown (the abscissa is the number of cycles Cycle, and the ordinate is the fluorescence value Fluorencence), if the PCR amplification reaction is successful and has a complete amplification curve (green solid line), then there are a baseline region, an exponential amplification region (the part circled by the purple circle), a linear growth region, a logarithmic growth region, and a saturation plateau region in this amplification curve. At the same time, during the actual amplification process, some non-specific signals or contaminations may be amplified in individual channels. At this time, the user can discard these abnormal data by setting a threshold by themselves. Therefore, there is also a manual threshold line (black dashed line) in the figure. Then, if an extension line (orange solid line) drawn from the linear segment of the amplification curve can be found, the intersection point of this extension line and the threshold line (blue five-pointed star in the figure) is behind the CT value determined by the manual threshold method (the value on the abscissa along the red vertical line of the pink diamond point, the point indicated by the black arrow). Considering that the amplification curve is also calculated by interpolation from discrete points, the method of using the extension line for calculation is less affected by the CT value calculation error caused by fluctuations in the exponential stage and is more accurate in principle.

[0043] However, the extension line can be accurately obtained on a continuous function but not in a discrete sequence. Calculating the intersection point with the threshold line in the exponential region has low calculation accuracy due to the small amount of collected data. Instead, it is better to use the fitting equation of the linear region. The intersection point of the fitting equation of the linear region and the baseline threshold line is more stable. The threshold line is obtained by parallelizing the threshold line based on the baseline fitting line at the starting segment of the amplification curve. In this way, the CT change caused by the baseline sinking is also removed. Solve for the intersection point using these two straight line equations, and use this intersection point as the CT value. This is the overall principle of obtaining the CT value in the present invention.

[0044] Based on the above principle, the present invention is elaborated in detail.

[0045] First, it is necessary to determine whether an effective amplification reaction has occurred. We collect the original fluorescence value amplification curve at the preset number of cycles during the quantitative polymerase chain reaction, fit the fluorescence values of the first preset number of cycles (cycles 1 to 12) into a straight line, and compare the difference with the corresponding original fluorescence amplification curve. When the difference between the fitted line and the original curve starts to continuously be greater than zero from a certain cycle point, it can be considered that an effective amplification reaction has occurred, and subsequent CT value calculation can be performed.

[0046] According to the principle proposed by the present invention, we need to obtain the ranges of the baseline region, the linear growth region, and the saturation plateau region, combined with Figure 2 It is not difficult to see that the exponential growth of the data between the baseline region and the linear growth region conforms to the characteristics of an elbow line. Therefore, there is an elbow point between the two regions, and the logarithmic growth of the data between the linear growth region and the saturation plateau region conforms to the characteristics of a knee line. Therefore, there is a knee point between the two regions. Then, as long as we find these two points, we can preliminarily determine the range of the linear region.

[0047] Based on this, we first define the following parameters using the original fluorescence values:

[0048]

[0049] In the formula, is the fluorescence value growth rate trend of the elbow point, is the fluorescence value change trend between cycle points, is the fluorescence value at the nth cycle. c is a ladder constant greater than or equal to 1. By examining different ladders, different curve characteristics will be obtained. When setting 2, it does not affect the sensitivity measurement for 45 cycle data. When setting 3, since the fluorescence value change range is small, the gap from the original curve characteristics increases. When setting 1, it is too sensitive to the fluctuations during the amplification process and is not conducive to extracting curve characteristics, but it has better sensitivity detection ability.

[0050] Since the numerical change trends of the baseline region, the linear growth region, and the saturation plateau region are relatively stable (stable without change or stable growth), the average fluorescence value between cycle points must be the fluorescence value at the midpoint between the two fluorescence values. Based on this, the closer the average value is to the fluorescence value at the midpoint, the closer it is to 0. Based on this, the preliminary typing of the amplification curve can be achieved through the search for points.

[0051] Based on the amplification curve after preliminary typing, we can screen for the elbow point in the exponential amplification region. Here, we utilize the characteristics that the fluorescence value of the elbow point is low and the change speed is fast, that is, as large as possible. By searching for points, the approximate position of the elbow point can be achieved acquisition

[0052] For the knee point, we first obtain the maximum fluorescence value M of the original fluorescence amplification curve and construct (where N is the system noise baseline), and define:

[0053]

[0054] In the formula, is the fluorescence value at the nth cycle after correction, is the growth rate trend of the fluorescence value at the knee point.

[0055] Then, according to the characteristics that the fluorescence value at the knee point is high and the change speed is fast, that is, as large as possible, by searching for the point position, the approximate position of the knee point is obtained.

[0056] Of course, considering whether the fitting is meaningful, we also consider the difference in the number of cycles between the approximate position of the elbow point and the approximate position of the knee point When the difference in the number of cycles is greater than 2, it is considered that the amplification reaction is normal and the CT value can be calculated. Otherwise, it proves that the curve has no obvious linear region and the amplification signal is weak or the peak appears too late. This difference is equivalent to measuring the difference between the baseline and the plateau of the amplification curve. If it is 2, the fitting loses meaning.

[0057] At this time, according to the approximate point of the elbow point and the approximate point of the knee point the linear region can be determined, and the fitting of the linear region equation and the baseline region equation can be carried out subsequently. Assume that the form of the fitted equation is y = kx + b, and discretize this straight-line equation:

[0058]

[0059] The original signal is denoted as , and thus the fitting error (i is a constant taking values from 1 to n) can be calculated, and the variance is calculated:

[0060]

[0061] When reaches the minimum value, there is an optimal solution. By solving this optimal solution, the fitting parameters and under the optimal fitting can be obtained.

[0062] For the fitting equation of the baseline region, the fitted equation is:

[0063]

[0064] Considering the threshold B set by the user, the final fitting equation for the baseline region is as follows:

[0065]

[0066] For the fitting equation of the linear region, the fitted equation is:

[0067]

[0068] In the formula, is the fitting parameter under the optimal fitting of the linear region equation, is the fitting parameter under the optimal fitting of the baseline region equation

[0069] The intersection point X of the two fitting equations can represent the elbow point position of the CT value.

[0070] Furthermore, in order to eliminate data jitter, after determining the linear region, data within the range of the number of cycles will also be smoothed. The data is divided by 100, rounded down, and then multiplied by 100.

[0071] In summary, for the automatic calculation method of the elbow point CT value of quantitative polymerase chain reaction described in the present invention, it segments using the characteristics of the line segment, performs fitting after segmentation, calculates the CT value using the characteristics of the fitted line, and its result fully reflects the characteristics of the original amplification curve. It has better compliance characteristics and more realistic effects for low-concentration samples such as weak positives. At the same time, the CV of the relative quantitative result is more stable.

[0072] It should be noted that all directional indications (such as up, down, left, right, front, back...) in the embodiments of the present invention are only used to explain the relative position relationship and movement conditions between components in a specific posture (as shown in the drawings). If the specific posture changes, the directional indications will also change accordingly.

[0073] In addition, in the present invention, descriptions such as "first", "second", and "one" are only for descriptive purposes and cannot be understood as indicating or implying their relative importance or implicitly indicating the quantity of the indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include at least one such feature. In the description of the present invention, "a plurality" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0074] In the present invention, unless otherwise clearly specified or defined, terms such as "connection" and "fixation" shall be understood in a broad sense. For example, "fixation" may be a fixed connection, a detachable connection, or integrated; it may be a mechanical connection or an electrical connection; it may be directly connected or indirectly connected through an intermediate medium, and it may be the communication inside two components or the interaction relationship between two components, unless otherwise clearly defined. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0075] In addition, the technical solutions between various embodiments of the present invention can be combined with each other, but it must be based on the ability of those of ordinary skill in the art to implement. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

Claims

1. A method for automatically calculating the elbow point CT value of quantitative polymerase chain reaction, characterized in that: Includes steps: S1: Obtain the original fluorescence value amplification curve at the preset cycle number during the quantitative polymerase chain reaction process; S2: Perform a preliminary judgment on the amplification curve according to the fluorescence value of the previous preset cycle, and proceed to the next step if it is determined that there is a valid amplification curve, otherwise end; S3: Perform amplification curve typing under linear judgment based on the original fluorescence value, and screen the approximate points of the elbow point based on the curve change characteristics of the exponential amplification period of the amplification curve, and screen the approximate points of the knee point based on the curve change characteristics of the saturation platform period of the amplification curve; S4: fitting the original fluorescence value amplification curve between the elbow approximation point and the knee approximation point, and the equation based on the best fit is the linear region equation; S5: fitting the amplification curve based on the original fluorescence value before the elbow point approximation point, and the equation based on the best fit is the baseline area equation; S6: Solve the intersection point based on the linear region equation and the baseline region equation, and use the intersection point as the target elbow point CT value; In the step S3, the amplification curve typing includes a baseline region, an exponential amplification region, a linear growth region, a logarithmic growth region and a saturation platform region; The amplification curve is typed by searching for the start and end points where the fluorescence value change trend between cycle points is stable. The fluorescence value change trend between cycle points is determined by the following formula: In the formula, is the fluorescence value variation trend between cycle points, c is the ladder constant greater than or equal to 1, is the fluorescence value at the nth cycle, by searching The corresponding cycle points are used to obtain the start and end points; The curve change characteristic of the exponential amplification period of the amplification curve is that the low fluorescence value state in the baseline area enters the exponential growth area where the fluorescence value exponentially increases. Among them, the fluorescence value growth rate trend of the elbow point is determined by the following formula: In the formula, The fluorescence value growth rate trend of the elbow point is obtained by searching the initial low fluorescence value in the original fluorescence value amplification curve and then the exponential growth range. The corresponding circulation points are used to screen the approximate points of the elbow points; The curve change characteristics of the amplification curve saturation platform period are that the fluorescence value in the linear area changes from linear growth to logarithmic growth area and finally enters the saturation platform area with stagnant growth. The growth rate trend of the fluorescence value at the knee point is determined by the following formula: Where M is the maximum fluorescence value, N is the system noise baseline, is the fluorescence value at the nth cycle after correction, The fluorescence value growth rate trend of the knee point is obtained by searching the fluorescence value logarithmically in the original fluorescence value amplification curve and finally entering the stagnant growth range. The corresponding cycle points are screened for knee point approximation points.

2. The method for automatically calculating the elbow point CT value of a quantitative polymerase chain reaction according to claim 1, characterized in that: In the step S2, the fluorescence values ​​of the preset cycles before fitting are formed into a straight line, and the difference is compared with the original fluorescence amplification curve corresponding to the part. When the difference is continuously greater than zero starting from any point, it is determined that there is a valid amplification curve.

3. The method for automatically calculating the elbow point CT value of a quantitative polymerase chain reaction according to claim 1, characterized in that: In the step S4, the linear region equation after fitting is expressed as: In the step S5, the fitted baseline region equation is expressed as: In the formula, x is the number of cycles, Y is the fluorescence value, is the fitting parameter of the optimal fitting of the linear region equation, is the fitting parameter of the optimal fitting of the baseline area equation, and B is the user-defined threshold.

Images