A method and device for analyzing thyroid cancer based on real-time fluorescent quantitative PCR

Abstract

This application discloses a method and apparatus for thyroid cancer analysis based on real-time quantitative PCR. The method includes: acquiring data to be analyzed; determining sample-specific threshold line data based on the data to be analyzed; calculating a corrected Ct value based on dual-inflection-point kinetics based on the sample-specific threshold line data to obtain the final corrected Ct value; and generating a detection result based on the final corrected Ct value. This application significantly improves the linearity of the fluorescence amplification curve, reduces the interference of background noise on the results, and improves the stability of Ct value calculation through optimized curve fitting algorithm processing.

Description

Technical Field

[0001] This application relates to the field of data processing technology, specifically to a method and apparatus for thyroid cancer analysis based on real-time quantitative PCR. Background Technology

[0002] The key performance indicator (Ct) for real-time quantitative PCR (polymerase chain reaction) technology is the cycle threshold (Ct), which refers to the number of PCR cycles required for the fluorescence signal in the reaction tube to first reach or exceed a set threshold. The accuracy of the Ct value directly affects the accuracy of PCR interpretation. A larger amount of initial template results in a faster fluorescence signal increase and a lower Ct value; conversely, a smaller amount of initial template results in a slower fluorescence signal increase and a higher Ct value.

[0003] The Thyroid Cancer 12-Gene Detection Kit (Fluorescent PCR Method) enables rapid qualitative detection of point mutations and fusion variants in thyroid cancer-related genes. This kit is used in conjunction with a fully automated medical PCR analysis system (SLAN-96S, National Medical Device Registration Certificate No. 20183221659). After setting the baseline and threshold according to the product instructions, the kit can interpret the positive and negative results for each locus.

[0004] However, in practical applications, the raw fluorescence amplification curves output by the SLAN-96S device often exhibit significant nonlinear distortion (manifested as poor curve smoothness, fluorescence signal drift, or fluctuation), leading to insufficient accuracy for technicians in setting baselines and thresholds. This can result in subjective errors and human interpretation differences, thereby increasing the risk of false negatives and false positives, potentially causing misjudgments in medical decisions and losses in patient management.

[0005] Currently, the conventional methods for calculating Ct values ​​in real-time quantitative PCR technology mainly include the threshold method and the maximum second derivative method.

[0006] 1. The threshold method is the most commonly used method for calculating Ct values. Its basic principle is to set a fixed fluorescence threshold line. When the fluorescence signal of the PCR amplification curve first exceeds this threshold, the corresponding cycle number is the Ct value. This method relies on the rapid rise of the fluorescence signal during the exponential amplification phase.

[0007] 2. The auxiliary baseline method calculates the average fluorescence value in the early cycles of the amplification curve (usually the first 3–15 cycles) as a baseline value. This baseline value is then subtracted from the entire amplification curve to correct for background signal, thereby improving the accuracy of Ct values.

[0008] 3. The maximum second derivative method calculates the second derivative of the amplification curve and finds the cycle number corresponding to its maximum value, which is then used as the Ct value. This method utilizes the curvature change of the amplification curve during the exponential growth phase, enabling a more objective determination of the Ct value and reducing the bias caused by manually setting thresholds.

[0009] The aforementioned traditional methods perform reasonably well under general experimental conditions, but they have significant shortcomings in the following scenarios:

[0010] 1. When the fluorescence signal is weak or the sample size is small, traditional methods are highly sensitive to the setting of thresholds and baselines, and are easily affected by noise, which leads to a significant decrease in the accuracy of Ct value calculation.

[0011] 2. When the amplification curve has obvious nonlinear characteristics (such as large fluctuations in the curve slope and high background noise), it is difficult to objectively determine the threshold using traditional methods, which leads to strong subjectivity in the interpretation by technicians and affects the stability and reliability of the test results.

[0012] 3. Furthermore, the nonlinearity of the amplification curve in the raw data output of SLAN-96S is particularly prominent, which exacerbates the limitations of traditional methods when applied to this system.

[0013] In summary, existing technical solutions are unable to effectively address the nonlinear distortion problem of PCR amplification curves, which severely restricts the objectivity and accuracy of test result interpretation and poses potential risks to actual medical diagnosis. Summary of the Invention

[0014] The purpose of this invention is to provide a method for analyzing thyroid cancer based on real-time quantitative PCR to at least solve one of the above-mentioned technical problems.

[0015] One aspect of the present invention provides a method for analyzing thyroid cancer based on real-time quantitative PCR, the method comprising:

[0016] Obtain the data to be analyzed;

[0017] Determine the sample-specific threshold line data based on the data to be analyzed;

[0018] Based on the sample-specific threshold line data, the correction Ct value is calculated based on the dual inflection point dynamics to obtain the final correction Ct value.

[0019] The detection results are generated based on the final corrected Ct value.

[0020] Optionally, the thyroid cancer analysis method based on real-time quantitative PCR includes:

[0021] The acquisition of the data to be analyzed includes:

[0022] Obtain the raw fluorescence signal from the PCR instrument for 1-35 cycles;

[0023] Based on the raw fluorescence signals of the positive control during 1-35 cycles of PCR, a smooth raw curve was fitted using the UnivariateSpline function.

[0024] Optionally, determining the sample-specific threshold data based on the data to be analyzed includes:

[0025] The difference between the maximum and minimum values ​​of the original fluorescence signal from 1 to 35 cycles is taken and multiplied by a coefficient of 0.05 as the calculation threshold. When the calculation threshold is less than 200, the calculation threshold is set to 200. Thus, the threshold on which the recalculation of the Ct value of the FAM / VIC channel corresponding to each reaction well site depends is obtained as the sample-specific threshold data.

[0026] Optionally, the step of calculating the correction Ct value based on the sample-specific threshold line data, thereby obtaining the final correction Ct value, includes:

[0027] Based on the original curve of each sample, the minimum value between 1 and 35 cycles is taken, and the sample-specific threshold line data is added to it as the signal inflection point. The cycle number corresponding to the inflection point is taken as the final correction Ct value.

[0028] Optionally, acquiring the data to be analyzed includes:

[0029] Obtain raw fluorescence data from 1-35 cycles of the PCR instrument;

[0030] Acquire instrument status logs, sample matrix parameters, and standard calibration data;

[0031] Extract features from the raw fluorescence data to obtain the raw fluorescence signal features;

[0032] Obtain the trained LSTM drift prediction model;

[0033] The features of the raw fluorescence data are input into the trained LSTM drift prediction model, which outputs a drift-free theoretical fluorescence signal.

[0034] The original fluorescence signal characteristics are corrected based on the drift-free theory of fluorescence signal to obtain the corrected fluorescence signal matrix;

[0035] Interference compensation is performed on the calibrated fluorescence signal matrix based on the sample matrix parameters to obtain the final fluorescence signal matrix as the data to be analyzed.

[0036] Optionally, the step of correcting the original fluorescence signal features based on the drift-free theory fluorescence signal to obtain the corrected fluorescence signal matrix includes:

[0037] The relative drift deviation is obtained from the fluorescence signal based on the drift-free theory.

[0038] A drift correction coefficient curve is generated based on the relative drift deviation;

[0039] A corrected fluorescence signal matrix is ​​generated based on the drift correction coefficient curve and the original fluorescence data.

[0040] Optionally, the relative drift deviation obtained from the fluorescence signal based on the drift-free theory is calculated using the following formula:

[0041] δ(x)= [F stdraw (x)-F stdpred (x)] / F stdpred (x);

[0042] Among them, F stdpred (x) The drift-free theoretical fluorescence signal of the x-th PCR reaction cycle; F stdraw (x) represents the raw fluorescence data of the x-th PCR reaction cycle; δ(x) represents the relative drift deviation of the x-th PCR reaction cycle.

[0043] Optionally, determining the sample-specific threshold data based on the data to be analyzed includes:

[0044] Extract positive control signals from the data to be analyzed;

[0045] Signal purification and feature extraction were performed on the positive control signal to obtain a noise-free ideal amplification curve;

[0046] Calculate the dynamic baseline threshold based on the noise-free ideal amplification curve;

[0047] Gene-specific threshold fine-tuning is performed on the dynamic baseline threshold to obtain gene-specific dynamic threshold lines as sample-specific threshold line data.

[0048] Optionally, the step of calculating the correction Ct value based on the sample-specific threshold line data, thereby obtaining the final correction Ct value, includes:

[0049] Intelligent baseline values ​​are generated based on the data to be analyzed and the sample-specific threshold line data.

[0050] The data to be analyzed is subjected to background subtraction using intelligent baseline values ​​to obtain the net signal sequence after background subtraction.

[0051] Generate a cycle-by-cycle amplification efficiency sequence based on the net signal sequence after background subtraction;

[0052] The fitting parameters are generated based on the cycle-by-cycle amplification efficiency sequence;

[0053] Based on the fitting parameters and the cycle-by-cycle amplification efficiency sequence, a Ct value compensation sequence for the critical cycle is generated.

[0054] Dual inflection point location information is obtained based on the net signal sequence after background subtraction.

[0055] An initial Ct value is generated based on the dual inflection point positioning information and the data to be analyzed.

[0056] The final corrected Ct value is obtained based on the initial Ct value and the Ct value compensation sequence of the key cycle.

[0057] This application also provides a thyroid cancer analysis device based on real-time quantitative PCR, the thyroid cancer analysis device based on real-time quantitative PCR comprising:

[0058] A data acquisition module for acquiring data to be analyzed;

[0059] A sample-specific threshold line data acquisition module is used to determine sample-specific threshold line data based on the data to be analyzed.

[0060] The final correction Ct value acquisition module is used to calculate the correction Ct value based on the sample-specific threshold line data, thereby obtaining the final correction Ct value.

[0061] The detection result acquisition module is used to generate detection results based on the final corrected Ct value.

[0062] The thyroid cancer analysis method based on real-time quantitative PCR proposed in this application has the following advantages:

[0063] 1. Objective and automated threshold and baseline determination: No longer relying on manually set thresholds and baselines, but automatically and dynamically adjusting according to the curve shape, eliminating the subjective bias of human interpretation;

[0064] 2. High-precision curve fitting optimization: Through the optimization of the curve fitting algorithm, the linearity of the fluorescence amplification curve is significantly improved, the interference of background noise on the results is reduced, and the stability of Ct value calculation is improved.

[0065] 3. High accuracy and stability verification: 1. Using the results of manual interpretation by experienced technicians as a reference, in the comparison of the results of the 12-gene detection kit for 1408 samples, the consistency between the manual interpretation results and the automated interpretation results of this application reached more than 98%.

[0066] 4. Reduced medical risks and improved interpretation efficiency: Significantly reduces the probability of false positives and false negatives, improving the safety and reliability of clinical diagnosis, while reducing reliance on operator experience and improving interpretation efficiency. The time per 96 wells is reduced from approximately 10-20 minutes to less than 2 minutes, greatly improving work efficiency. Attached Figure Description

[0067] Figure 1 This is a schematic flowchart of a thyroid cancer analysis method based on real-time quantitative PCR according to an embodiment of this application.

[0068] Figure 2 This is a comparative schematic diagram of the present application and SLAN, which is an embodiment of the present application.

[0069] Figure 3 This is a schematic diagram of the automatic analysis results of an embodiment of this application. Detailed Implementation

[0070] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions in the embodiments of this application will be described in more detail below with reference to the accompanying drawings. In the drawings, the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The described embodiments are some, but not all, embodiments of this application. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain this application, and should not be construed as limiting this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application. The embodiments of this application will be described in detail below with reference to the accompanying drawings.

[0071] Example 1:

[0072] like Figure 1 The thyroid cancer analysis method based on real-time quantitative PCR shown includes:

[0073] Step 1: Obtain the data to be analyzed;

[0074] Step 2: Determine the sample-specific threshold line data based on the data to be analyzed;

[0075] Step 3: Calculate the correction Ct value based on the sample-specific threshold line data to obtain the final correction Ct value;

[0076] Step 4: Generate the detection results based on the final corrected Ct value.

[0077] In this embodiment, acquiring the data to be analyzed includes:

[0078] Obtain the raw fluorescence signal from 1-35 cycles of the PCR instrument; specifically, first extract the raw fluorescence signal from the positive control from 1-35 cycles of stage 4 of the PCR instrument.

[0079] Based on the raw fluorescence signals of the positive control during 1-35 cycles of PCR, a smooth raw curve was fitted using the UnivariateSpline function.

[0080] In this embodiment, determining the sample-specific threshold data based on the data to be analyzed includes:

[0081] For each site, the difference between the maximum and minimum values ​​of the corrected fluorescence signal between cycles 1 and 35 is taken and multiplied by a factor of 0.05 to obtain the calculated threshold line (Δcut-off). When Δcut-off < 200, Δcut-off is set to 200. Thus, the threshold line on which the recalculated Ct value of the FAM / VIC channel corresponding to each reaction well site depends is obtained as the max and min values ​​of the corrected fluorescence signal for the positive control within cycles 1-35.

[0082] In this embodiment, the step of calculating the correction Ct value based on the sample-specific threshold line data, thereby obtaining the final correction Ct value, includes:

[0083] Based on the FAM / VIC channel threshold lines determined for the corresponding reaction wells in the positive control, the corrected Ct value is obtained by adding the minimum Y-axis value of the original curve to the X-axis value (cycle number) corresponding to the threshold line. Specifically, a smooth original curve is fitted using the UnivariateSpline function based on the original fluorescence signal of each sample in PCR instrument stage 4 from 1 to 35 cycles. The minimum fluorescence signal between 1 and 35 cycles is taken, and the threshold line is added to it as the signal inflection point. The cycle number corresponding to the inflection point is taken as the corrected Ct value.

[0084] Step 3: Interpretation of test results: Based on the corrected Ct value, interpret the results using the positive judgment value specified in the instruction manual.

[0085] The thyroid cancer analysis method based on real-time quantitative PCR in this application may also include the following:

[0086] Plot the amplification curve obtained after correcting the original fluorescence signal, and label the corrected Ct value and the corresponding corrected fluorescence intensity for easy manual verification. Specifically, set the baseline start and end cycle numbers to 1-10, subtract the baseline from the original signal, and then use a four-parameter logic model (4PL model, see Formula 1) and the least squares method to correct the original fluorescence signal. Simultaneously, the original curve, the corrected amplification curve, and the amplification curve obtained from the PCR instrument are presented for easy manual verification.

[0087] Formula 1:

[0088] Where a is the maximum asymptote (plateau value), d is the minimum asymptote (baseline value), c is the x value at the inflection point (maximum slope point), b is the Hill slope, which determines the steepness of the curve's inflection, and x is the cycle number.

[0089] Meanwhile, the client system adopts the Python Dash framework and B / S architecture to implement the following technical solutions:

[0090] Interactive data upload and processing engine. Supports drag-and-drop upload of raw PCR data in Excel / CSV format, dynamically calls backend algorithms (including the four core steps mentioned above), automatically calculates Ct values ​​and generates mutation detection results, and supports multi-threaded batch processing.

[0091] Real-time visualization system. Utilizing the Dash DataTable component, it dynamically renders data, supports adaptive column widths, multi-line text folding, and scrolling, ensuring clear display of the results.

[0092] Auditable data flow management. A dual-folder system isolates raw data and results (uploads / results), and file names embed timestamps to ensure data traceability and version control. An encrypted download interface ( / download / ) is provided via Flask routing. <filename>This meets the requirements for clinical data compliance.

[0093] See Figure 3 The method described in this application allows for the rapid and automated interpretation of negative, positive, or manually verified results. The analysis results include the project name, sample name, interpretation result (positive, negative, PC abnormal, NTC normal, retest), and detailed information (site, reaction well location, Ct value). A download button is also provided to download detailed results (including corrected amplification curves).

[0094] See Figure 2 Because SLAN uses a fixed threshold of 200, the Ct value at the TP53 site in this sample was 26.8, resulting in a positive interpretation. However, the method in this application uses a unique threshold for each mutation site, derived from the difference between the maximum and minimum values ​​of the calibrated fluorescence signal curve of the positive control sample at that site, multiplied by a coefficient of 0.05 (200 is used when the threshold is <200). In this case, the threshold for the TP53 mutation site obtained using this method was 323.4. Therefore, the Ct value at the TP53 site in this sample exceeded 35, resulting in a preliminary negative result. Since the difference in Ct values ​​between the two samples was >2, this situation could be annotated and manually verified. After retesting, the TP53 site result was negative, and the patient's postoperative outcome was also benign. Furthermore, no significant amplification was observed from the original curve. However, because SLAN uses a fixed threshold, it could lead to incorrect interpretation and potentially misdiagnosis.

[0095] Example 2 (only the parts that differ from Example 1 are described):

[0096] In Example 2, acquiring the data to be analyzed includes:

[0097] Obtain raw fluorescence data from 1-35 cycles of the PCR instrument; in this embodiment, the raw fluorescence data are specifically as follows: the dimension is the number of reaction wells × 35 cycles, including the test sample, positive / negative control, and 3 concentration standards (10). 3 / 10 4 / 10 5 FAM / VIC channel fluorescence values ​​(copies / μL)

[0098] Acquire instrument status logs (e.g., heating module temperature T(x) per cycle, excitation light path intensity L(x), photomultiplier tube gain (G(x))), sample matrix parameters (e.g., nucleic acid purity of the sample to be tested (A260 / A280 nucleic acid purity ratio, representing the absorbance of the sample at 260 nm (nucleic acid absorption peak) and 280 nm (protein absorption peak) as detected by a UV spectrophotometer, and the ratio of the two is calculated), and inhibitor concentration, representing the concentration of substances in the sample that can inhibit Taq DNA polymerase activity (such as hemoglobin in blood, humic acid in tissue samples, and residual phenolic reagents during extraction) (C). inh (ng / μL) and standard calibration data (theoretical Ct values ​​of 3 standards);

[0099] Extract features from the raw fluorescence data to obtain the raw fluorescence signal features;

[0100] Obtain the trained LSTM drift prediction model;

[0101] The features of the raw fluorescence data are input into the trained LSTM drift prediction model, which outputs a drift-free theoretical fluorescence signal.

[0102] The original fluorescence signal characteristics are corrected based on the drift-free theory of fluorescence signal to obtain the corrected fluorescence signal matrix;

[0103] Interference compensation is performed on the calibrated fluorescence signal matrix based on the sample matrix parameters to obtain the final fluorescence signal matrix as the data to be analyzed.

[0104] In this embodiment, the trained LSTM drift prediction model can be trained in the following way:

[0105] Using historical datasets as the training set (70%) and validation set (30%), mean squared error (MSE) as the loss function and Adam as the optimizer, the model was iteratively trained until the validation set MSE < 0.001, and the trained model was saved.

[0106] In this embodiment, the historical dataset can be obtained in the following way:

[0107] A dataset was created using raw fluorescence signals of standard samples from 100 historical batches of experiments, along with instrument status parameters (300 sets in total, 35 cycles per set).

[0108] The theoretical C values ​​of three standards in historical data (used to back-calculate the ideal fluorescence signal label).

[0109] Based on the theoretical Ct value of the standard, and according to the PCR amplification kinetic model (exponential growth + plateau phase), the ideal fluorescence signal (signal value under drift-free state) for each cycle is deduced.

[0110] In this embodiment, the specific structure of the LSTM model is as follows:

[0111] Input layer: 4D (original fluorescence signal F of standard) stdraw (x), T (x), L (x), G (x));

[0112] Hidden layers: 2 layers, 128 neurons per layer, ReLU activation function;

[0113] Output layer: 1D (ideal fluorescence signal F) stdideal (x));

[0114] In this embodiment, the original fluorescence signal F of the three standards from the current batch (the batch to be analyzed) is used. stdraw Inputting (x), T(x), L(x), and G(x)) into a trained LSTM model in cyclic order, and outputting the drift-free theoretical fluorescence signal F for each cycle. stdpred (x);

[0115] In this embodiment, the step of correcting the original fluorescence signal features based on the drift-free theory fluorescence signal to obtain the corrected fluorescence signal matrix includes:

[0116] The relative drift deviation is obtained based on the fluorescence signal according to the drift-free theory; in this embodiment, the relative drift deviation is calculated using the following formula:

[0117] δ(x)= [F stdraw (x)-F stdpred (x)] / F stdpred (x);

[0118] Among them, F stdpred (x) The drift-free theoretical fluorescence signal of the x-th PCR reaction cycle; F stdraw δ(x) represents the raw fluorescence data of the x-th PCR reaction cycle; δ(x) represents the relative drift deviation of the x-th PCR reaction cycle; (δ(x)>0 indicates signal overestimation, δ(x)<0 indicates signal underestimation).

[0119] A drift correction coefficient curve is generated based on the relative drift deviation; specifically, the CubicSpline function from the Python scipy library is used to smoothly interpolate δ(x) over 35 loops to generate a continuous drift correction coefficient curve K. drift (x).

[0120] A corrected fluorescence signal matrix is ​​generated based on the drift correction coefficient curve and the original fluorescence data; specifically, the original fluorescence value F of each reaction well is calculated according to the cycle number (x=1-35). raw (x) and the corresponding cycle K drift (x) Match;

[0121] For each cycle, the fluorescence value of all reaction wells was adjusted for drift using the formula: F driftcorr (x)= F raw (x) / K drift (x);

[0122] The corrected signals were organized in a cycle format of reaction well number × 35 to form a corrected fluorescence signal matrix.

[0123] In this embodiment, interference compensation is performed on the calibration fluorescence signal matrix based on the sample matrix parameters to obtain the final fluorescence signal matrix as the data to be analyzed, including:

[0124] Obtain the matrix interference coefficient weights; in this embodiment, the matrix interference coefficient weights can be obtained in the following way:

[0125] The coefficients obtained by fitting a large amount of preliminary experimental data (500+ clinical samples) represent A260 / A280 and C. inh Contribution to interference with PCR signal;

[0126] Specifically, through multiple linear regression analysis, different A260 / A280 ratios and different C values ​​were analyzed. inh Samples of different concentrations were subjected to PCR amplification, and their signal attenuation rates were recorded. The values ​​were then calculated using a regression model, which in this embodiment were 0.12 and 0.08.

[0127] For each sample to be tested, the matrix interference coefficient K is calculated according to the formula. matrix The specific formula is as follows:

[0128] K matrix =1+0.12×(2.0-A260 / A280)+0.08×lg(C inh +1).

[0129] Based on the correspondence between the sample and the reaction well, the F of each reaction well was... driftcorr (x) and the K of the corresponding samples matrix match;

[0130] The drift-corrected signal is compensated cycle by cycle to obtain the final fluorescence signal for each reaction well, as shown in the following formula:

[0131] F finalcorr (x)=F driftcorr (x)×K matrix ;

[0132] The individual fluorescence signals are combined to form a final fluorescence signal matrix, which serves as the data to be analyzed.

[0133] In this embodiment, determining the sample-specific threshold data based on the data to be analyzed includes:

[0134] Extract positive control signals from the data to be analyzed;

[0135] In this embodiment, extracting the positive control signal from the data to be analyzed includes:

[0136] Based on the well location markings in the experimental metadata, locate the row corresponding to the positive control in the two-dimensional corrected fluorescence signal matrix (e.g., mark well B1 as the FAM channel positive control, then find the row corresponding to well B1 in the matrix).

[0137] From the located rows, the corrected fluorescence values ​​of 35 cycles (x=1 to 35) were extracted to form the positive control corrected signal sequence F. ctrlcorr (x), whose data format is [F ctrlcorr (1), F ctrlcorr (2), ...,F ctrlcorr (35)];

[0138] If the detection involves both FAM and VIC channels, repeat the above two steps: extract the positive control correction signal sequence for each channel based on the positive control well markings for both channels, i.e., F... ctrlcorr FAM(x) (FAM channel) and F ctrlcorrVIC (x) (VIC channel).

[0139] Signal purification and feature extraction were performed on the positive control signal to obtain a noise-free ideal amplification curve;

[0140] In this embodiment, signal purification and feature extraction are performed on the positive control signal to obtain a noise-free ideal amplification curve, including:

[0141] The positive control correction signal sequence F ctrlcorr (x) Input the Python pywt library and perform a 3-level db4 wavelet decomposition;

[0142] After decomposition, four components were obtained: low-frequency approximation component A3 (representing the core trend of positive control amplification, without noise), high-frequency detail component D1 (corresponding to electronic noise in cycles 1-5, showing high-frequency random fluctuations), high-frequency detail component D2 (corresponding to reagent background noise in cycles 6-15, showing mid-frequency stable fluctuations), and high-frequency detail component D3 (corresponding to signal attenuation noise in cycles 16-35, showing a slow decrease in low frequency).

[0143] D1 (Electronic Noise) Suppression: Kalman filtering is used, with the state equation set as A3(x) = A3(x-1) + w(x) (where w(x) is the process noise with a variance of 0.02), and the observation equation as F ctrlcorr (x) = A3(x) + D1(x); Substituting A3 and D1 into the filtering model, the signal after eliminating electronic noise is calculated and denoted as A3 + D2 + D3. filtered .

[0144] D2 (reagent background) inhibition: Based on the stage 4 parameters, cycles 6-15 were determined as the main reagent background interval. A3+D2+D3 within this interval were calculated. filtered The mean value of the signal is denoted as B (i.e., the reagent background baseline value); a background exponential decay compensation function B is constructed. comp (x) = B × exp (-0.05 × (x - 15)) (x ≥ 6, because the reagent background decays exponentially with the amplification process after 15 cycles); A3 + D2 + D3 filtered Subtract B comp (x), the signal after eliminating the reagent background is obtained, denoted as A3+D3. compensated .

[0145] D3 (signal attenuation) suppression: Calculate A3+D3 within 16-35 cycles. compensated Signal attenuation coefficient α=[ln(F ctrlcorr (35) / F ctrlcorr (16))] / (35-16) (reflects the signal attenuation rate); construct the attenuation compensation function D3 comp (x)=F ctrlcorr (x)×exp(-α×(x-16)); put A3+D3 compensated Multiply by D3 comp (x), and finally obtain the noise-free ideal amplification curve, denoted as F. ctrlideal (x).

[0146] For noise-free ideal amplification curve F ctrlideal (x) Find the first derivative to obtain the amplification rate curve F'(x) = F ctrlideal (x)-F ctrlideal (x-1) (x≥2);

[0147] Find the cyclic interval where F'(x) ≥ 0.1 × max(F'(x)) (i.e., the stage where the amplification rate reaches more than 10% of the maximum rate), and denote it as the exponential growth segment, with its initial cycle being x. start The loop terminates at x. end (Usually a 10→25 cycle);

[0148] Extract the maximum value of F'(x) within this interval, and denote it as the maximum amplification rate S. max .

[0149] Calculate the dynamic baseline threshold based on the noise-free ideal amplification curve;

[0150] In this embodiment, calculating the dynamic baseline threshold based on the noise-free ideal amplification curve includes:

[0151] Calculate the base threshold Δcut-off base :

[0152] From the noise-free ideal amplification curve F ctrlideal In (x), the maximum fluorescence value (denoted as Max fluorescence) and the minimum fluorescence value (denoted as Min fluorescence) of cycles 1-35 are extracted.

[0153] Introducing an amplification rate weighting factor W(x): due to the exponential growth segment (x) start →x end The region ) is the core region for specific signal amplification. To improve threshold specificity, x is set at x... start →x end When W(x) is 1.2, W(x) = 1.0 for the rest of the cycles;

[0154] Calculate the base threshold using the following formula: Δcut-off base =(Ma fluorescence - Min fluorescence) × 0.05 × W(x).

[0155] Calculate the instrument stability coefficient K stab :

[0156] Based on the real-time temperature log T(x), the instrument stability coefficient for each cycle is calculated using the following formula: K stab (x)=1-|T(x)-T target | / T target ;

[0157] This coefficient reflects the degree of deviation between the instrument's actual temperature and the target temperature: K stab (x)=1 indicates no deviation (instrument stability), K stab The smaller (x) is, the greater the deviation (instrument instability), for example, T(x) = 60.6℃, T target At 60℃, K stab (x) = 1 - |60.6 - 60| / 60 = 0.99.

[0158] Calculate the final dynamic threshold Δcut-off final :

[0159] The base threshold is combined with the instrument stability coefficient using the following formula: Δcut-off temp = Δcut-off base ×K stab (x);

[0160] The final dynamic threshold for each cycle is Δcut-off. final (x)=max (Δcut-off temp ,200) (200 is the lower limit of the threshold);

[0161] For example, a certain loop Δcut-off temp =190 (below 200), then Δcut-off final (x)=200; if Δcut-off temp =220 (higher than 200), then Δcut-off final (x)=220.

[0162] Gene-specific threshold fine-tuning is performed on the dynamic baseline threshold to obtain gene-specific dynamic threshold lines as sample-specific threshold line data.

[0163] In this embodiment, gene-specific threshold fine-tuning is performed on the dynamic baseline threshold to obtain gene-specific dynamic threshold lines as sample-specific threshold line data, including:

[0164] Based on the name of the gene being detected (e.g., BRAF V600E), retrieve the ideal amplification rate range (e.g., 800-1000 fluorescence values / cycle) and ideal threshold fluctuation range (e.g., ±5%) of the corresponding gene from the 12-gene dynamics database (existing database).

[0165] Calculate the actual maximum amplification rate S of the current positive control. max实际 The ideal S of this gene in the database max理想均值 The deviation rate is calculated using the formula: Deviation rate = |S max实际 -S max理想均值 | / S max理想均值 ×100%.

[0166] If the deviation rate is ≤10% (the actual amplification rate is basically consistent with the ideal rate), then Δcut-off will not be adjusted. final (x), directly use this threshold;

[0167] If the deviation rate > 10%: If S max实际 >S max理想均值 (If amplification is too rapid, it may lead to false positives), then Δcut-off will be implemented. final (x) Increase by 3% (within the ideal threshold fluctuation range); if S max实际 max理想均值 (If amplification is too slow, it may lead to false negatives), then Δcut-off will be implemented. final ​(x) is reduced by 3% to ensure that the threshold is adapted to the current gene amplification characteristics.

[0168] In this embodiment, the step of calculating the correction Ct value based on the sample-specific threshold line data, thereby obtaining the final correction Ct value, includes:

[0169] Intelligent baseline values ​​are generated based on the data to be analyzed and the sample-specific threshold line data.

[0170] In this embodiment, generating a smart baseline value based on the data to be analyzed and the sample-specific threshold line data includes:

[0171] Locate the row of reaction wells corresponding to the current sample in the data to be analyzed;

[0172] The 35 cyclic fluorescence values ​​from this row were extracted to form the corrected fluorescence signal sequence F for a single test sample. samplecorr (x) (x = 1 - 35, x is the repeating number).

[0173] Preset filtering parameters (window size = 5, polynomial order = 2).

[0174] The Savitzky-Golay filtering algorithm is used, with 5 consecutive cycles as the sliding window (e.g., when calculating x=10, the window covers x=8-12).

[0175] For each window, F samplecorr The (x) value is fitted with a second-order polynomial, and the intermediate value of the fitted curve is used to replace the original signal value to obtain the filtered signal.

[0176] Iterate through loops 1-35 to generate the complete filtered signal sequence F. samplefilter (x).

[0177] Calculate the second derivative F'' of the loop 1-15 sample (x) (x≥3):

[0178] F'' sample (x)=[F samplefilter (x)-F samplefilter [(x-1)]-F samplefilter (x-1)-F samplefilter [(x-2)];

[0179] Iterate through loops 1-15 to find F's for 5 consecutive loops. ample (x) The interval where the absolute value of each x is less than or equal to the threshold (e.g., x = 2 - 6 in a loop);

[0180] Calculate F within this interval samplefilter The mean of (x) is used as the intelligent baseline value B. smart .

[0181] The data to be analyzed is subjected to background subtraction using intelligent baseline values ​​to obtain the net signal sequence after background subtraction.

[0182] In this embodiment, the process of performing background subtraction on the data to be analyzed using an intelligent baseline value to obtain a net signal sequence after background subtraction includes:

[0183] Background subtraction is performed on the filtered signal for each loop:

[0184] F samplesub (x)=F samplefilter (x)-B smart ;

[0185] If the calculation result is negative (low-cycle signal is below baseline), then F samplesub (x) is set to 0 (to avoid negative signals interfering with subsequent efficiency calculations).

[0186] Generate a cycle-by-cycle amplification efficiency sequence based on the net signal sequence after background subtraction;

[0187] In this embodiment, generating the cycle-by-cycle amplification efficiency sequence based on the net signal sequence after background subtraction includes:

[0188] Starting with the cycle x=2 (requiring the signal from the previous cycle), calculate the amplification efficiency cycle by cycle:

[0189] E(x) = [F samplesub (x)-F samplesub [(x-1)] / F samplesub (x-1)×100%;

[0190] If F samplesub (x-1)=0 (no net signal in the previous cycle), then E(x) is set to 0 (marked as amplification not started).

[0191] The fitting parameters are generated based on the cycle-by-cycle amplification efficiency sequence;

[0192] In this embodiment, generating fitting parameters based on the cycle-by-cycle amplification efficiency sequence includes:

[0193] An efficiency decay model is constructed using the Logistic function:

[0194] E pred (x) =E max / [1+exp(k×(x-x0))];

[0195] Where k is the decay rate coefficient and x0 is the efficiency decay start cycle, both of which are parameters to be fitted;

[0196] Using E(x) (x=2-35) as the measured value, k and x0 are optimized using the nonlinear least squares method to make E pred The mean square error (MSE) of E(x) and E(x) is less than 5%, generating a fitted efficiency decay curve sequence E pred (x).

[0197] Based on the fitting parameters and the cycle-by-cycle amplification efficiency sequence, a Ct value compensation sequence for the critical cycle is generated.

[0198] In this embodiment, generating the Ct value compensation sequence for the critical cycle based on the fitting parameters and the cycle-by-cycle amplification efficiency sequence includes:

[0199] Efficiency deviation calculated cycle by cycle: ΔE(x) = E max -E pred (x) (reflects the gap between actual efficiency and ideal efficiency);

[0200] Constructing a Ct value compensation model: Ct comptemp (x) = x - [ΔE(x) / 100] × 0.8;

[0201] Centered on x0 (the decay start cycle), extract Ct for x0±3 cycles. comptemp (x) (covering the possible Ct value ranges afterward), thereby obtaining the Ct value compensation sequence for the key cycle.

[0202] Dual inflection point location information is obtained based on the net signal sequence after background subtraction.

[0203] In this embodiment, the dual-inflection point localization information based on the net signal sequence after background subtraction includes:

[0204] Calculate F samplesub The second derivative F'' of (x) sub (x) (x≥3):

[0205] F'' sub (x)=[F samplesub (x)-F samplesub (x-1)]-[F samplesub (x-1)-F samplesub [(x-2)];

[0206] Positioning the acceleration inflection point Ct accel : Finding the first one to satisfy F'' sub (x)>0 and F'' sub The cycle of (x-1) ≤ 0 serves as the acceleration inflection point Ct. accel (The expansion transitioned from slow to exponential growth).

[0207] Locating the deceleration inflection point Ct_decel: Finding the first condition that satisfies F'' sub (x) < 0 and F'' sub The cycle of (x-1) ≥ 0 serves as the deceleration inflection point Ct_decel (the expansion transitions from exponential growth to a plateau).

[0208] If no inflection point is found, it is marked as no effective amplification detected.

[0209] An initial Ct value is generated based on the dual inflection point positioning information and the data to be analyzed.

[0210] In this embodiment, generating the initial Ct value based on the dual inflection point positioning information and the data to be analyzed includes:

[0211] Extracting Ct from gene-specific dynamic threshold lines accel Thresh (Ct) of the cycle accel ), Ct decel Thresh(Ct) threshold of the loop decel );

[0212] Validity judgment: If F samplesub (Ct accel )≥Thresh(Ct accel (If the acceleration inflection point signal has crossed the threshold, ensuring effective amplification), then the initial Ct value is calculated using the following formula:

[0213] (Geometric mean equilibrium inflection point fluctuation);

[0214] If the signal does not cross the threshold, it is marked as the amplified signal not reaching the threshold, and the Ct value is invalid.

[0215] The final corrected Ct value is obtained based on the initial Ct value and the Ct value compensation sequence of the key cycle.

[0216] In this embodiment, obtaining the final corrected Ct value based on the initial Ct value and the Ct value compensation sequence of the key cycle includes:

[0217] Take Ct init The integer part of (initial Ct value) is used as the compensation loop x. comp (For example, Ct) init =23.92, then x comp =23);

[0218] From the cyclically changing temporary compensation Ct value sequence Ct comptemp In (x), the compensation loop x is invoked. comp The corresponding temporary compensation value Ct is denoted as Ct. comptemp (x comp For example, if x comp (Compensation loop = 23, call Ct) comptemp (23) (The temporary compensation loop threshold corresponding to loop 23 is 22.88;

[0219] The final corrected Ct value is calculated using the following formula:

[0220] Ct final =Ct init +(Ct comptemp (x comp )-x comp );

[0221] Among them, (Ct) comptemp (x comp )-x comp (This is to compensate for the deviation, such as 22.88-23=-0.12, which means compensating for -0.12 cycles).

[0222] If Ct final ≤38, output Ct final If Ct final If the value is >38 or Ct is invalid, it is marked as not detected.

[0223] This application has the following advantages:

[0224] To address the common issues of drift and matrix interference in raw fluorescence data from PCR instruments, a hierarchical correction system was constructed to improve signal quality at its source.

[0225] By using an LSTM drift prediction model and combining it with instrument status logs (temperature, optical path intensity, photomultiplier tube gain), the theoretically drift-free fluorescence signal is inferred. Then, the deviation per cycle is calculated using a relative drift deviation formula, generating a smooth drift correction coefficient curve. Compared to the traditional method of ignoring drift with a fixed threshold, this method can dynamically compensate for signal overestimation / underestimation issues in different cycles.

[0226] By introducing sample matrix parameters such as nucleic acid purity (A260 / A280) and inhibitor concentration (Cinh), and using weighting coefficients (0.12, 0.08) fitted from over 500 sets of previous clinical data, a matrix interference coefficient was constructed. This coefficient can effectively compensate for signal attenuation and avoid false negatives caused by differences in sample quality, particularly in scenarios involving hemoglobin interference in blood samples and humic acid residue in tissue samples.

[0227] A combined approach employing 3-layer db4 wavelet decomposition, Kalman filtering, and exponential decay compensation is used to accurately separate electronic noise (D1), reagent background noise (D2), and signal attenuation noise (D3) from the positive control signal, ultimately generating a noise-free ideal amplification curve. Compared to the simple noise reduction of traditional baseline subtraction, this approach preserves signal specificity while reducing noise interference.

[0228] By combining the maximum / minimum fluorescence values ​​of the noise-free ideal amplification curve, the weighting factor of the exponential growth segment, and the stability coefficient calculated from the instrument's real-time temperature deviation, a dynamic baseline threshold for each cycle is generated. For example, when the instrument temperature deviation is large, the threshold will be automatically adjusted to avoid misjudgment and adapt to the instrument status under different operating conditions.

[0229] Using a 12-gene kinetics database, the actual amplification rate (S) is calculated based on the ideal amplification rate range of the currently detected gene (e.g., BRAF V600E). max实际 ) and ideal speed (S) max 理想 The deviation rate is calculated. If the deviation rate is >10%, the threshold is adjusted up / down accordingly (e.g., 3% up for excessively fast amplification to prevent false positives, and 3% down for excessively slow amplification to prevent false negatives), to solve the problem of incompatibility of thresholds caused by differences in amplification efficiency of different genes.

[0230] The acceleration inflection point (Ct) is located by calculating the second derivative of the net signal after background subtraction. accel (Exponential expansion phase) and deceleration inflection point (Ct) decel (When amplification enters the plateau phase), the initial Ct value is calculated only when the acceleration inflection point signal crosses the threshold, avoiding misjudging non-specific amplification as a valid signal and reducing false positives.

[0231] Based on the logistic efficiency decay model, the cycle-by-cycle amplification efficiency is fitted to generate a Ct value compensation sequence for key cycles. Finally, by combining the initial Ct value with the compensation deviation, the Ct value shift caused by amplification efficiency fluctuations is corrected. In scenarios with low template amount samples (Ct close to 38), the calculation error can be controlled within ±0.2 cycles.

[0232] This application can effectively adapt to various scenarios, including samples with low nucleic acid purity (A260 / A280 deviation from 2.0), samples containing inhibitors (such as residual phenolic reagents), and scenarios with fluctuating instrument status (temperature deviation, optical path attenuation), by using stratified calibration to avoid detection failure due to sample or instrument problems.

[0233] For the application scenarios of the 12-gene detection kit for thyroid cancer, by fine-tuning the gene-specific threshold, it can simultaneously process the amplification curves of different genes (such as BRAF, RAS, RET), without having to set parameters for each gene separately, thus improving detection efficiency and ensuring high consistency of results.

[0234] This application also provides a thyroid cancer analysis device based on real-time quantitative PCR, which includes a data acquisition module, a sample-specific threshold line data acquisition module, a final correction Ct value acquisition module, and a detection result acquisition module; wherein,

[0235] The data acquisition module is used to acquire the data to be analyzed.

[0236] The sample-specific threshold line data acquisition module is used to determine the sample-specific threshold line data based on the data to be analyzed.

[0237] The final correction Ct value acquisition module is used to calculate the correction Ct value based on the sample-specific threshold line data, thereby obtaining the final correction Ct value.

[0238] The detection result acquisition module is used to generate detection results based on the final corrected Ct value.

[0239] Although the present invention has been described in detail above with general descriptions and specific embodiments, modifications or improvements can be made to it, which will be obvious to those skilled in the art. Therefore, all such modifications or improvements made without departing from the spirit of the present invention fall within the scope of protection claimed by the present invention.< / filename>

Claims

1. A real-time fluorescent quantitative PCR-based thyroid cancer analysis device, characterized by, The thyroid cancer analysis device based on real-time quantitative PCR includes: The data acquisition module is used to acquire data to be analyzed; the acquisition of data to be analyzed includes: Obtain raw fluorescence data from 1-35 cycles of the PCR instrument; Acquire instrument status logs, sample matrix parameters, and standard calibration data; Extract features from the raw fluorescence data to obtain the raw fluorescence signal features; Obtain the trained LSTM drift prediction model; The features of the raw fluorescence data are input into the trained LSTM drift prediction model, which outputs a drift-free theoretical fluorescence signal. The original fluorescence signal characteristics are corrected based on the drift-free theory of fluorescence signal to obtain the corrected fluorescence signal matrix; Interference compensation is performed on the calibrated fluorescence signal matrix based on the sample matrix parameters to obtain the final fluorescence signal matrix as the data to be analyzed; The step of correcting the original fluorescence signal features based on the drift-free theory fluorescence signal to obtain the corrected fluorescence signal matrix includes: The relative drift deviation is obtained from the fluorescence signal based on the drift-free theory. A drift correction coefficient curve is generated based on the relative drift deviation; A corrected fluorescence signal matrix is ​​generated based on the drift correction coefficient curve and the original fluorescence data; The relative drift deviation obtained from the fluorescence signal based on the drift-free theory is calculated using the following formula: δ(x) = [F stdraw (x) - F stdpred (x)] / F stdpred (x) ; where F stdpred (x) is the non-drifted theoretical fluorescence signal for the xth PCR reaction cycle; F stdraw (x) is the raw fluorescence data for the xth PCR reaction cycle; δ(x) is the relative drift bias for the xth PCR reaction cycle. A sample-specific threshold data acquisition module, wherein the sample-specific threshold data acquisition module is used to determine sample-specific threshold data based on the data to be analyzed; the determination of sample-specific threshold data based on the data to be analyzed includes: Extract positive control signals from the data to be analyzed; Signal purification and feature extraction were performed on the positive control signal to obtain a noise-free ideal amplification curve; Calculate the dynamic baseline threshold based on the noise-free ideal amplification curve; Gene-specific threshold fine-tuning is performed on the dynamic baseline threshold to obtain gene-specific dynamic threshold lines as sample-specific threshold line data; The final correction Ct value acquisition module is used to calculate the correction Ct value based on dual-inflection-point dynamics according to the sample-specific threshold line data, thereby obtaining the final correction Ct value; the calculation of the correction Ct value based on dual-inflection-point dynamics according to the sample-specific threshold line data to obtain the final correction Ct value includes: Intelligent baseline values ​​are generated based on the data to be analyzed and the sample-specific threshold line data. The data to be analyzed is subjected to background subtraction using intelligent baseline values ​​to obtain the net signal sequence after background subtraction. Generate a cycle-by-cycle amplification efficiency sequence based on the net signal sequence after background subtraction; The fitting parameters are generated based on the cycle-by-cycle amplification efficiency sequence; Based on the fitting parameters and the cycle-by-cycle amplification efficiency sequence, a Ct value compensation sequence for the critical cycle is generated. Dual inflection point location information is obtained based on the net signal sequence after background subtraction. An initial Ct value is generated based on the dual inflection point positioning information and the data to be analyzed. The final corrected Ct value is obtained based on the initial Ct value and the Ct value compensation sequence of the key cycle; The detection result acquisition module is used to generate detection results based on the final corrected Ct value.

Images