Correlation analysis method for gas chromatography detection data

Abstract

The invention relates to the field of data processing, in particular to a correlation analysis method for gas chromatography detection data, which comprises the following steps of: extracting peak shape characteristics of each peak, and constructing peak shape characteristic penalty terms among different peaks according to peak shape characteristic differences; and calculating the DTW distance of the local chromatograms to which the corresponding data points belong at the corresponding time point pairs by combining the difference condition between the corresponding signal intensities of the local chromatograms to which the different data points are located at the time point pairs and a dynamic time warping algorithm. Aiming at the peak overlapping problem of a standard DTW algorithm in gas chromatographic data correlation analysis, the sensitivity of the DTW algorithm to the peak shape characteristics is improved by introducing a self-adaptive multi-characteristic fusion penalty term based on the peak shape characteristics in distance calculation, so that the contribution of different compounds in overlapped peaks can be distinguished more accurately, and the accuracy of the gas chromatographic data correlation analysis is improved. Therefore, the correlation between the compounds is evaluated more accurately, and the influence of peak overlapping on a correlation analysis result is reduced.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of data processing, and particularly relates to a correlation analysis method for gas chromatography detection data. BACKGROUND

[0002] Gas chromatography technology is widely used in environmental monitoring, food safety, drug analysis, petrochemical industry and other fields, and is an important means for qualitative and quantitative analysis of substances. In the analysis of complex samples, gas chromatography data usually contains peaks of multiple compounds, and there may be correlations between these compounds. For example, some compounds may be intermediates or products on the same reaction path, or they have the same source in the sample or are affected by similar environmental factors. Analyzing the correlation between these compounds helps to reveal the source of substances in the sample, transformation process, interaction, etc., which is of great significance for sample component analysis, pollution tracing, metabolic pathway research, etc.

[0003] At present, the commonly used methods for analyzing the correlation of different compounds in gas chromatography data mainly include: directly calculating the Pearson correlation coefficient, Spearman correlation coefficient, etc. of the peak area or peak height of different compounds, or using cross-correlation analysis method to calculate the cross-correlation coefficient of two chromatograms (or local chromatograms). In recent years, dynamic time warping (DTW) algorithm has also been introduced into chromatography data analysis for comparing the similarity of different chromatograms or local chromatograms. The DTW algorithm can better handle the shape difference of peaks and the retention time drift problem by allowing nonlinear alignment of the time axis, thereby improving the accuracy of correlation analysis to some extent.

[0004] Although the above methods have been applied to a certain extent in the correlation analysis of gas chromatography data, they all have a significant limitation: in the presence of peak overlap, the accuracy of correlation calculation will be severely affected. In gas chromatography data, due to the complexity of sample components, instrument condition fluctuations, etc., the peaks of different compounds often partially or completely overlap in retention time. When the peaks overlap, the directly measured peak area, peak height or local chromatogram signal intensity is actually the superposition of multiple compound signals, and cannot accurately reflect the true content or trend of a single compound. This leads to unreliable correlation analysis results based on these measured values, which may mask the true correlation between compounds or produce false correlations.

[0005] Although the standard DTW algorithm can handle a certain degree of peak shape change and retention time drift, the difference in signal intensity is mainly considered when calculating the distance between two local chromatograms, and the peak shape feature information is not fully utilized. In the peak overlapping area, the difference in signal intensity is actually a mixture of the differences in signals of multiple compounds, and the standard DTW algorithm cannot distinguish whether this difference is caused by the true difference between compounds or by peak overlap. This leads to distortion of the distance calculated by the DTW algorithm in the case of peak overlap, and further affects the accuracy of the correlation analysis. SUMMARY

[0006] The purpose of the present application is to extract a correlation analysis method for gas chromatography detection data to solve the problems raised in the background art.

[0007] In order to solve the above technical problems, the present application extracts the following technical scheme: a correlation analysis method for gas chromatography detection data, the method comprising:

[0008] S1, acquiring signal intensity values arranged in time sequence and retention times corresponding to each signal intensity value by gas chromatography data acquisition by a gas chromatograph, constructing original data points, and preprocessing the original data points; performing peak recognition on the preprocessed data using a watershed algorithm, dividing the chromatogram into multiple independent peak regions, and extracting local chromatograms corresponding to each peak;

[0009] S2, extracting peak shape features of each peak, constructing peak shape feature penalty terms between different peaks according to the difference in peak shape features, and calculating the DTW distance between the local chromatograms corresponding to the corresponding data points at the corresponding time points based on the difference between the signal intensities at the corresponding time points and the dynamic time warping algorithm;

[0010] S3, dynamically adjusting the weights of different features in the peak shape feature penalty term according to the characteristics of the local chromatogram to which each data point belongs, optimizing the DTW distance calculation formula of the local chromatogram corresponding to the corresponding data points at the corresponding time points based on the optimized weights of different features in the peak shape feature penalty term, and calculating the DTW distance between all local chromatograms of all peaks; converting the obtained DTW distance into a correlation coefficient to construct a correlation coefficient matrix;

[0011] S4, based on the constructed correlation coefficient matrix, dynamically screening the obtained peak pairs, and determining the correlation between the compounds represented by the obtained peak pairs according to the screening results.

[0012] The dynamic time warping (DTW) algorithm is a commonly used time series similarity measurement method in gas chromatography data correlation analysis. The DTW algorithm can better handle the shape difference of peaks and the retention time drift problem by allowing nonlinear alignment of the time axis. To some extent, the accuracy of correlation analysis is improved. The standard DTW algorithm only considers the difference in signal intensity when calculating the distance between two local chromatograms, ignoring the peak shape feature information. In the presence of peak overlap, the signal intensity of the local chromatogram cannot accurately reflect the characteristics of a single compound because the signals of different compounds are mixed together. If the DTW distance is calculated only according to the signal intensity difference, it will lead to incorrect alignment results, which will further affect the accuracy of the correlation analysis. For example, two local chromatograms may have a large difference in signal intensity due to overlapping peaks of different compounds, but if these overlapping peaks are decomposed, it will be found that the peak shape characteristics of some compounds are very similar, while the peak shape characteristics of other compounds are quite different. The standard DTW algorithm cannot distinguish this situation and considers that the two local chromatograms are not similar, thereby underestimating the correlation between certain compounds.

[0013] Further, in the process of extracting the local chromatogram corresponding to each peak in S1, the vertex coordinates of each peak are extracted, and the first element of the vertex coordinates is the corresponding retention time, and the second element is the peak height; according to the peak identification result, taking the vertex coordinates of each peak as the center, each peak is extracted to form the local chromatogram corresponding to the peak.

[0014] Further, the peak shape features in S2 include peak change trend, peak sharpness and peak shape skewness.

[0015] In the process of constructing the peak shape feature penalty term between different peaks in S2, the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinates at time point t1 is denoted as s i (t1);The signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinates at time point t2 is denoted as s j (t2);

[0016] The peak change trend corresponding to s i (t1) is denoted as s i ‘(t1);The peak change trend corresponding to s j (t2) is denoted as s j ‘(t2);

[0017] The peak sharpness corresponding to s i (t1) is denoted as s i ‘’(t1);The peak sharpness corresponding to s j (t2) is denoted as s j “(t2);

[0018] s i (t1) the corresponding peak shape skewness is denoted as sw i (t1);s j (t2) the corresponding peak shape skewness is denoted as sw j (t2);

[0019] s i (t1) and s j (t2) between the peak shape feature penalty term is denoted as P(i, t1, j, t2);

[0020] P(i,t1,j,t2)=exp(a·Δf (i,t1,j,t2) )-1

[0021] Wherein, a represents a preset adjustment factor; Δf (i,t1,j,t2) ) represents the integrated peak shape feature difference between s i (t1) and s j (t2);

[0022]

[0023] Wherein, The n th peak shape feature weight calculated by the signal intensity of the local chromatogram corresponding to the peak pair with data point i as the vertex coordinate at time point t1, and each The initial value of each of s

[0024] When n=1, And

[0025] When n=2, And

[0026] When n=3, And

[0027] The application introduces a penalty term based on peak shape features in the DTW distance calculation; Specifically, when calculating the distance between two local chromatograms, not only the difference in signal intensity is considered, but also the difference in peak shape features; If the peak shape feature difference of two local chromatograms on a certain time point pair is large, then the distance between them is increased to reduce the possibility of aligning these two points by the DTW algorithm. In this way, the DTW algorithm can be more inclined to align points with similar peak shape features, thereby improving its alignment accuracy in the case of peak overlap, and further improving the accuracy of correlation analysis.

[0028] Further, the s i (t1) and s j (t2) are obtained in the following wayi (t1) and s j The first derivative of (t2) is as follows:

[0029] s′ i (t1)=(s i (t1+1)-s i (t1-1)) / 2

[0030] s′ j (t2)=(s j (t2+1)-s j (t2-1)) / 2

[0031] Among them, s i (t1+1) represents the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1+1; s i (t1-1) represents the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1-1; s j (t2+1) represents the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2+1; s j (t2-1) represents the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2-1;

[0032] The s i ''(t1) and s j The method of obtaining (t2) (s) i (t1) and s j The second derivative of (t2) is as follows:

[0033] s″ i (t1)=(s i (t1+1)-2*s i (t1)+s i (t1-1))

[0034] s″ j (t2)=(s j (t2+1)-2*s j (t2)+s j (t2-1))

[0035] The sw i (t1) and sw j (t2) is the quotient obtained by dividing the third-order center distance corresponding to the signal intensity of the data point in the corresponding local chromatogram by the cube of the standard deviation corresponding to the signal intensity of the data point in the corresponding local chromatogram.

[0036] Further, the S2 calculates the DTW distance of the local chromatogram to which the corresponding data point belongs at the corresponding time point pair, involving the following calculation formula:

[0037] d ij (t1,t2)=|s i (t1)-s j (t2)|+λ·P(i,t1,j,t2)

[0038] Wherein, d ij (t1,t2) represents the DTW distance between s i (t1) and s j (t2); |s i (t1)-s j (t2)| represents the signal intensity difference of the two local chromatograms corresponding to s i (t1) and s j (t2) at the corresponding time points; λ represents the preset penalty term weight coefficient for adjusting the intensity of the peak shape feature penalty term.

[0039] Further, the S3 dynamically adjusts the weight of different features in the peak shape feature penalty term according to the characteristics of the local chromatogram to which each data point belongs, involving the following calculation formula:

[0040]

[0041] In the process of optimizing the DTW distance calculation formula of the local chromatogram to which the corresponding data point belongs at the corresponding time point pair based on the optimized weight of different features in the peak shape feature penalty term, each is obtained when n is different respectively, and then substituted into the calculation formula of the comprehensive peak shape feature difference between s i (t1) and s j (t2), to obtain the updated result of the DTW distance between s i (t1) and s j (t2).

[0042] The weight of different peak shape features in the penalty term in step S2 of the present application is fixed, i.e., a pre-set constant. The problem of this fixed weight mode is that it is assumed that the importance of all peak shape features for distinguishing different compounds is the same, and this importance is consistent at all data points and all time points. However, in actual gas chromatography data, the importance of different peak shape features can vary with data points and time. For example, in the top region of the peak, the signal intensity is large, and the signal-to-noise ratio is high, at this time the peak height and the peak shape skewness can better reflect the characteristics of the compound, and the importance of the first derivative and the second derivative (which mainly reflect the edge characteristics of the peak) is relatively low. While in the edge region of the peak, the signal intensity is small, and the signal-to-noise ratio is low, at this time the reliability of the peak height and the peak shape skewness decreases, and the importance of the first derivative and the second derivative (which reflect the rising and falling trend of the peak) becomes relatively high. In addition, the importance of different peak shape features can also be different for different types of compounds.

[0043] In order to solve the problem caused by the fixed weight, the present application constructs an adaptive dynamic weight mechanism in this step, i.e., dynamically adjusts the weight of different features in the peak shape feature penalty term according to the characteristics of the local chromatogram of each data point. Specifically, the weight of each peak shape feature is calculated according to the signal intensity, the absolute value of the first derivative, and the absolute value of the second derivative of each data point. Through this adaptive dynamic weight mechanism, the peak shape feature penalty term is more targeted and flexible, further improving the correlation analysis accuracy of the DTW algorithm in the case of peak overlap.

[0044] Further, the calculation formula for converting the obtained DTW distance into a correlation coefficient in S3 is as follows:

[0045] r ij = exp(-d DTW / σ)

[0046] wherein d DTW represents the DTW distance corresponding to any one peak pair formed by the local chromatograms of all peaks; σ represents an adjustment parameter, and the value of σ is equal to the average value of the DTW distances corresponding to each peak pair; r ij represents the correlation coefficient corresponding to d DTW .

[0047] The correlation coefficient matrix is the summary result of the correlation coefficients corresponding to each peak pair, and the correlation coefficient in the g1th row and the g2th column in the correlation coefficient matrix represents the correlation coefficient between the g1th peak and the g2th peak, wherein the value range of g1 and g2 is [1, H], and H represents the total number of identified peaks.

[0048] Further, S4 includes:

[0049] Calculate the mean μ and the standard deviation σ1 of all non-diagonal elements in the correlation coefficient matrix, and calculate the screening threshold μ+k*σ for dynamic screening of the obtained peak pairs, wherein k represents a preset adjustable parameter;

[0050] Identify the peak pairs corresponding to the elements greater than the screening threshold in the correlation coefficient matrix as the relevant peak pairs, and determine that there is a correlation between the compounds represented by the relevant peak pairs; identify the peak pairs corresponding to the elements less than or equal to the screening threshold in the correlation coefficient matrix as the irrelevant peak pairs, and determine that there is no correlation between the compounds represented by the irrelevant peak pairs.

[0051] Compared with the prior art, the beneficial effects achieved by the present application are: the present application proposes an improved DTW distance calculation method aiming at the peak overlap problem existing in the correlation analysis of gas chromatography data based on the standard DTW algorithm, introduces an adaptive multi-feature fusion penalty term based on peak shape features in distance calculation, improves the sensitivity of the DTW algorithm to peak shape features, enables it to more accurately distinguish the contributions of different compounds in overlapping peaks, and thus more accurately evaluates the correlation between compounds and reduces the influence of peak overlap on the correlation analysis results. BRIEF DESCRIPTION OF DRAWINGS

[0052] The accompanying drawings are used to extract a further understanding of the present application, and constitute a part of the specification, together with the embodiments of the present application, for explaining the present application, and do not constitute a limitation on the present application. In the drawings:

[0053] Figure 1 is a flowchart of a correlation analysis method for gas chromatography detection data according to the present application. DETAILED DESCRIPTION

[0054] The technical solutions in the embodiments of the present application will be described in detail below with reference to the accompanying drawings of the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0055] Please refer to Figure 1 The technical solution of the present application is: a correlation analysis method for gas chromatography detection data, the method comprises:

[0056] S1, gas chromatography data acquisition is performed by a gas chromatograph, signal intensity values arranged in time sequence and retention times corresponding to each signal intensity value are obtained, to form original data points, and the original data points are preprocessed; the preprocessed data is subjected to a watershed algorithm for peak identification, the chromatogram is divided into a plurality of independent peak regions, and the local chromatogram corresponding to each peak is extracted;

[0057] In the process of extracting the local chromatogram corresponding to each peak in S1, the vertex coordinates of each peak are extracted, the first element of the vertex coordinates being the corresponding retention time and the second element being the peak height; according to the peak identification result, a local chromatogram corresponding to each peak is constructed by taking the vertex coordinates of each peak as the center and taking m data points before and after the vertex coordinates.

[0058] The preprocessing process in this embodiment includes: smoothing filtering of the original data by using a Savitzky-Golay filter to remove high-frequency noise in the data, the filter window size and the polynomial order being determined according to the actual data characteristics; then, baseline correction is performed by using an iterative weighted least squares method to eliminate baseline drift caused by changes in the state of the instrument or the presence of interfering substances in the sample, the number of iterations and the weight function being determined according to the actual data characteristics; then, abnormal data points are detected and removed by setting a threshold and a local outlier factor (LOF).

[0059] In S2, the peak shape features of each peak are extracted, and a peak shape feature penalty term between different peaks is constructed according to the differences in the peak shape features; the DTW distance of the local chromatogram to which the corresponding data point belongs at the corresponding time point pair is calculated by combining the differences between the signal intensities at the time point pair of the local chromatogram where the corresponding data point is located and the dynamic time warping algorithm.

[0060] The peak shape features in S2 include peak change trend, peak sharpness and peak shape skewness.

[0061] In the process of constructing the peak shape feature penalty term between different peaks in S2, the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinates at time point t1 is denoted as s i (t1), t1 being the retention time of data point i; the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinates at time point t2 is denoted as s j (t2), t2 being the retention time of data point j;

[0062] The peak change trend corresponding to s i (t1) is denoted as s i ‘(t1); the peak change trend corresponding to s j (t2) is denoted as s j ‘(t2); the acquisition methods of s i ‘(t1) and s j ‘(t2) are as follows:

[0063] s′ i (t1) = (s i (t1+1) - s i(t1-1)) / 2

[0064] s′ j (t2) = (s j (t2+1) - s j (t2-1)) / 2

[0065] wherein s i (t1+1) represents the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1+1; s i (t1-1) represents the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1-1; s j (t2+1) represents the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2+1; s j (t2-1) represents the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2-1;

[0066] s i (t1) is recorded as s i '(t1); s j (t2) is recorded as s j '(t2); the s i '(t1) and s j '(t2) are obtained in the following manner:

[0067] s″ i (t1) = (s i (t1+1) - 2*s i (t1) + s i (t1-1))

[0068] s″ j (t2) = (s j (t2+1) - 2*s j (t2) + s j (t2-1))

[0069] s i (t1) is recorded as sw i (t1); s j (t2) is recorded as sw j (t2);

[0070] The sw i (t1) and sw j (t2) are each obtained by dividing the third central moment corresponding to the signal intensity of the data points in the corresponding local chromatogram by the cube of the standard deviation corresponding to the signal intensity of the data points in the corresponding local chromatogram.

[0071] The peak shape characteristic penalty term between s i (t1) and s j (t2) is denoted as P(i, t1, j, t2), which is used to measure the difference in peak shape characteristics between two local chromatograms, and when the difference in peak shape characteristics is large, the value of P(i, t1, j, t2) is large, thereby increasing the distance between the two local chromatograms and reducing the possibility of matching them when they are aligned by DTW;

[0072] P(i, t1, j, t2) = exp(a·Δf (i,t1,j,t2) )-1

[0073] wherein a represents a preset adjustment factor, and in the embodiment, the value range of a is generally [0.5, 2], and the specific value is determined by cross-validation; the greater the value of a, the more significant the influence of the difference in peak shape characteristics on the penalty term; Δf (i,t1,j,t2) ) represents the comprehensive difference in peak shape characteristics between s i (t1) and s j (t2), which comprehensively considers the difference in multiple peak shape characteristics;

[0074]

[0075] wherein, represents the nth peak shape characteristic weight calculated by the signal intensity of the local chromatogram corresponding to the peak pair with data point i as the vertex coordinate at time point t1, and each of the different values of n corresponds to a different peak shape characteristic; The initial value of each of the different values of n is 1.

[0076] When n = 1, and

[0077] When n = 2, and

[0078] When n = 3, and

[0079] The S2 calculates the DTW distance of the local chromatogram to which the corresponding data point belongs at the corresponding time point pair, and the calculation formula involved is as follows:

[0080] d ij (t1, t2) = |s i (t1) - s j (t2)| + λ·P(i, t1, j, t2)

[0081] wherein d ij (t1, t2) represents the distance between s i (t1) and sj (t2) between the two, which is composed of two parts: signal intensity difference (Euclidean distance) and peak shape feature penalty term; |s i (t1)-s j (t2)| represents the signal intensity difference (Euclidean distance) of the two local chromatograms corresponding to s i (t1) and s j (t2) respectively at the corresponding time point; λ represents the preset penalty term weight coefficient for adjusting the intensity of the peak shape feature penalty term; the value range of λ is [0.1, 1], and the specific value is determined by cross-validation. The greater the value of λ, the greater the influence of the peak shape feature difference on the final distance; the smaller the value of λ, the smaller the influence of the peak shape feature difference.

[0082] In this embodiment, it is assumed that there are two local chromatograms A and B corresponding to two data points in the gas chromatography data. Both of the two local chromatograms contain a chromatographic peak, but there is a certain difference in peak shape and there is partial overlap. Specifically, it is assumed that the peak of A is sharp and the retention time is early, and the peak of B is flat and the retention time is late, and there is partial overlap between the falling edge of the peak of A and the rising edge of the peak of B.

[0083] If only the traditional Euclidean distance is used to calculate the DTW distance, the DTW algorithm will tend to align the peak top of A with the rising edge or a point in the overlapping region of the peak of B, because in this way the sum of the signal intensity difference can be minimized. This alignment method is obviously unreasonable because it ignores the significant difference in shape of the two peaks. However, the method proposed in the present application introduces a peak shape feature penalty term when calculating the DTW distance. The penalty term considers the first derivative difference, the second derivative difference and the peak shape skewness difference of the two local chromatograms at the corresponding time point. Since the peak of A is sharp, the absolute value of its second derivative is large; the peak of B is flat, and the absolute value of its second derivative is small. Therefore, when the DTW algorithm tries to align the peak top of A with the rising edge or the overlapping region of the peak of B, the second derivative difference will increase significantly, resulting in an increase in the value of the peak shape feature penalty term and further increasing the final distance.

[0084] Similarly, due to the asymmetry of the peak shapes of A and B, the peak shape skewness difference will contribute to the penalty term. Considering the influence of these peak shape feature differences, the DTW algorithm will tend to align the peak top of A with the peak top of B, because although this alignment method will cause the signal intensity difference to increase slightly, it can significantly reduce the peak shape feature difference, so that the final distance containing the penalty term reaches a minimum value.

[0085] S3, dynamically adjusting the weights of different features in the peak shape feature penalty term according to the features of the local chromatogram to which each data point belongs, optimizing the DTW distance calculation formula of the local chromatogram to which the corresponding data point belongs at the corresponding time point pair based on the optimized weights of different features in the peak shape feature penalty term, and calculating the corresponding DTW distance between each pair of peaks; converting the obtained DTW distance into a correlation coefficient to construct a correlation coefficient matrix;

[0086] The S3 dynamically adjusts the weights of different features in the peak shape feature penalty term according to the features of the local chromatogram to which each data point belongs, and the calculation formula is as follows:

[0087]

[0088] In this embodiment, an adaptive dynamic weight mechanism is introduced to dynamically adjust the weights according to the features of the local chromatogram of each data point. Specifically, the greater the signal intensity, the greater the weight. This is because the signal at the peak is more reliable and can better represent the characteristics of the peak. The greater the local change rate (absolute value of first derivative), the smaller the weight, because the absolute value of the first derivative is large at the rising and falling edges of the peak, but these regions have less contribution to the peak shape similarity; the greater the sharpness of the peak (absolute value of second derivative), the smaller the weight, because the absolute value of the second derivative is large near the peak vertex, but these regions have less contribution to the peak shape similarity; in this way, the dynamic weight can automatically adjust the contribution of different peak shape features in the penalty term according to the features of the local chromatogram.

[0089] To clearly explain the principle and effect of the dynamic weight mechanism proposed in this step, the assumptions used in step are still used. Near the peak top: the signal intensity of A and B is large, but the absolute value of the second derivative of A is large (the peak is sharp), and the absolute value of the second derivative of B is small (the peak is flat); according to the dynamic weight calculation formula, the second peak shape feature weight of A will be smaller (because the absolute value of the second derivative of A is large), and the second peak shape feature weight of B will be larger, which means that the penalty term pays more attention to the difference between the peak height and the peak shape skewness near the peak top, reducing the influence of the difference in the second derivative (the difference in the sharpness of the peak), because the peak height can better reflect the characteristics of the compound near the peak top, while the sharpness of the peak is easily affected by noise and other factors.

[0090] Edge of the peak: the signal intensity of A and B is small, but the first derivative absolute value of A is large (the peak is steep), and the first derivative absolute value of B is small (the peak is flat). According to the dynamic weight calculation formula, the first peak shape feature weight of A will be small (because the first derivative absolute value of A is large), and the first peak shape feature weight of B will be large, which means that at the edge of the peak, the penalty term pays more attention to the difference in peak shape skewness, and reduces the influence of the difference in first derivative (the difference in peak steepness), because the signal intensity is small at the edge of the peak, and the signal-to-noise ratio is low, resulting in a decrease in the reliability of the first derivative, while the peak shape skewness can better reflect the overall shape of the peak; in this way, the dynamic weight can adaptively adjust the contribution of different peak shape features in the penalty term according to the characteristics of the local chromatogram, so that the penalty term is more targeted and flexible. Thus, the alignment accuracy of the DTW algorithm under the condition of peak overlap is further improved, and finally the accuracy of the correlation analysis is improved.

[0091] In the process of optimizing the DTW distance calculation formula of the local chromatogram corresponding to the data points belonging to the corresponding time point pair based on the optimized weights of different features in the peak shape feature penalty term, the obtained respective is substituted into the calculation formula of the comprehensive peak shape feature difference between s i (t1) and s j (t2), to obtain the update result of the DTW distance between s i (t1) and s j (t2). The calculation formula for converting the obtained DTW distance into a correlation coefficient in S3 is as follows:

[0092] r ij = exp(-d DTW / σ)

[0093] where d DTW represents the DTW distance of any peak pair formed by the local chromatograms of all peaks; σ represents an adjustment parameter, and the value of σ is equal to the average of the DTW distances of the respective peak pairs; r ij represents the correlation coefficient corresponding to d DTW .

[0094] The correlation coefficient matrix constructed is the summary result of the correlation coefficients of the obtained peak pairs, and the correlation coefficient between the g1th peak and the g2th peak in the g1th row and the g2th column of the correlation coefficient matrix represents the correlation coefficient between the g1th peak and the g2th peak, and the value range of g1 and g2 is [1, H], and H represents the total number of identified peaks.

[0095] S4, based on the constructed correlation coefficient matrix, dynamically screening the obtained peak pairs, and determining the correlation between the compounds represented by the obtained peak pairs according to the screening result.

[0096] The S4 comprises:

[0097] The mean μ and the standard deviation σ1 of all non-diagonal elements in the correlation coefficient matrix are calculated, and a screening threshold μ+k*σ for dynamically screening the obtained peak pairs is calculated, wherein k represents a preset adjustable parameter, and is usually 2 or 3, the greater the value of k, the fewer the screened correlation peak pairs but the stronger the correlation; the smaller the value of k, the more the screened correlation peak pairs but some weakly correlated peak pairs are included; the specific value of k should be adjusted in combination with the actual sample complexity and analysis requirements;

[0098] The peak pairs corresponding to the elements greater than the screening threshold in the correlation coefficient matrix are identified as the correlation peak pairs, and it is determined that there is correlation between the compounds represented by the correlation peak pairs; the peak pairs corresponding to the elements less than or equal to the screening threshold in the correlation coefficient matrix are identified as the non-correlation peak pairs, and it is determined that there is no correlation between the compounds represented by the non-correlation peak pairs.

[0099] It should be noted that, in this document, the terms such as first and second are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "contain" or any other variants thereof are intended to cover non-exclusive inclusion, so that the process, method, article or equipment including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such process, method, article or equipment.

[0100] Finally, it should be noted that: the above only describes the preferred embodiments of the present application, and does not limit the present application, although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement to some technical features. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method for correlation analysis of gas chromatography detection data, characterized by, The method comprises: S1, acquiring signal intensity values arranged in time sequence and corresponding retention times of each signal intensity value by gas chromatography data acquisition through a gas chromatograph, to constitute original data points, and preprocessing the original data points; adopting a watershed algorithm to perform peak identification on the preprocessed data, segmenting a chromatogram into multiple independent peak regions, and extracting a local chromatogram corresponding to each peak; S2, extracting peak shape features of each peak, constructing a peak shape feature penalty term between different peaks according to the peak shape feature differences, combining the difference between the corresponding signal intensities at each time point pair of the local chromatogram where the data points are located and the dynamic time warping algorithm to calculate the DTW distance of the local chromatogram where the corresponding data points belong at the corresponding time point pair; S3, dynamically adjusting the weights of different features in the peak shape feature penalty term according to the features of the local chromatogram where each data point belongs, optimizing the DTW distance calculation formula of the local chromatogram where the corresponding data points belong at the corresponding time point pair based on the weights of different features in the optimized peak shape feature penalty term, and calculating the DTW distance between all peak local chromatograms to form each peak pair; converting the obtained DTW distance into a correlation coefficient to construct a correlation coefficient matrix; S4, based on the constructed correlation coefficient matrix, dynamically screening the obtained peak pairs, and determining the correlation between the compounds represented by the obtained peak pairs according to the screening results.

2. The method for correlation analysis of gas chromatography detection data according to claim 1, characterized in that: In the process of extracting the local chromatogram corresponding to each peak in S1, the vertex coordinates of each peak are extracted, and the first element in the vertex coordinates is the corresponding retention time and the second element is the peak height; according to the peak identification result, taking the vertex coordinates of each peak as the center, each peak corresponds to a local chromatogram composed of m data points in front and behind.

3. The method for correlation analysis of gas chromatography detection data according to claim 1, characterized in that: The peak shape features in S2 include peak change trend, peak sharpness and peak shape skewness; In the process of constructing the peak shape feature penalty term between different peaks in S2, the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1 is denoted as s i (t1); the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2 is denoted as s j (t2). s i (t1) the corresponding peak change trend is denoted as s i '(t1); s j (t2) the corresponding peak change trend is denoted as s j '(t2); s i (t1) the corresponding peak sharpness is noted s i (t1); s j (t2) the corresponding peak sharpness is noted s j (t2); s i (t1) the corresponding peak shape skewness is noted sw i (t1); s j (t2) the corresponding peak shape skewness is noted sw j (t2); s i (t1) and s j (t2) is denoted by P(i, t1, j, t2); P(i, t1, j, t2) = exp(a - Af (i,t1,j,t2) )-1 wherein a represents a preset adjustment factor; Δf (i,t1,j,t2) ) represents the difference between the integrated peak shape characteristics of s i (t1) and s j (t2). wherein, represents the nth peak shape feature weight calculated from the signal intensity of the local chromatogram corresponding to the peak with vertex coordinate of data point i at time point t1, and n corresponds to each of the different values The initial value of each of the peak shape feature weights is 1. when n = 1, and when n = 2, and when n = 3, and 4. The method for correlation analysis of gas chromatography detection data according to claim 3, characterized in that: The s i ‘(t1) and s j ‘(t2) are obtained as follows: s′ i (t1) = (s i (t1+1)-s i (t1-1)) / 2 s′ j (t2) = (s j (t2+1)-s j (t2-1)) / 2 where s i (t1+1) represents the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1+1; s i (t1-1) represents the signal intensity of the local chromatogram corresponding to the peak with data point i as the vertex coordinate at time point t1-1; s j (t2+1) represents the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2+1; s j (t2-1) represents the signal intensity of the local chromatogram corresponding to the peak with data point j as the vertex coordinate at time point t2-1; The s i The s j The s s" i (t1) = (s i (t1+1)-2*s i (t1)+s i (t1-1)) s" j (t2) = (s j (t2+1)-2*s j (t2)+s j (t2-1)) The sw i (t1) is the quotient of the third central moment divided by the cube of the standard deviation of the signal intensities of the data points within the respective local chromatogram. j (t2) is the quotient of the third central moment divided by the cube of the standard deviation of the signal intensities of the data points within the respective local chromatogram.

5. The method for correlation analysis of gas chromatography detection data according to claim 3, wherein: The calculation formula involved in calculating the DTW distance of the local chromatogram where the corresponding data points belong at the corresponding time point pair in S2 is as follows: d ij (t1,t2) = |s i (t1)-s j (t2)| + λ · P(i, t1, j, t2) wherein d ij (t1, t2) represents the DTW distance between s i (t1) and s j (t2); |s i (t1) - s j (t2) | represents the signal intensity difference of the two local chromatograms corresponding to s i (t1) and s j (t2) at the corresponding time points, respectively; and λ represents a preset penalty term weight coefficient for adjusting the intensity of the peak shape feature penalty term.

6. The method for correlation analysis of gas chromatography detection data according to claim 3, wherein: In S3, the calculation formula involved in dynamically adjusting the weights of different features in the peak shape feature penalty term according to the features of the local chromatogram where each data point belongs is as follows: In the process of optimizing the DTW distance calculation formula of the local chromatogram pair to which the corresponding data points belong at the corresponding time point pair based on the weights of different features in the optimized peak shape feature penalty term, the respective is substituted into the calculation formula of the comprehensive peak shape feature difference between s i (t1) and s j (t2), and the update result of the DTW distance between s i (t1) and s j (t2) is obtained.

7. The method for correlation analysis of gas chromatography detection data according to claim 1, wherein: The calculation formula for converting the obtained DTW distance into a correlation coefficient in S3 is as follows: r ij = exp(_d DTW / σ) wherein d DTW represents the DTW distance between any two peaks of the local chromatogram; σ represents an adjustment parameter, and the value of σ is equal to the average of the DTW distances corresponding to each peak pair; r ij represents d DTW the corresponding correlation coefficient; The correlation coefficient matrix constructed is the summary result of the correlation coefficients of each peak pair, and the element in the g1th row and the g2th column in the correlation coefficient matrix represents the correlation coefficient between the g1th peak and the g2th peak, and the value range of g1 and g2 is [1, H], and H represents the total number of identified peaks.

8. The method for correlation analysis of gas chromatography detection data according to claim 1, wherein: S4 comprises: Calculate the mean μ and standard deviation σ1 of all non-diagonal elements in the correlation coefficient matrix, and calculate the screening threshold μ+k*σ for dynamically screening the obtained peak pairs, wherein k represents a preset adjustable parameter; The peak pair corresponding to the element greater than the screening threshold in the correlation coefficient matrix is identified as a correlation peak pair, and it is determined that there is correlation between the compounds represented by the correlation peak pair; the peak pair corresponding to the element less than or equal to the screening threshold in the correlation coefficient matrix is identified as an unrelated peak pair, and it is determined that there is no correlation between the compounds represented by the unrelated peak pair.

Images