Method and system for processing methylation sequencing data

Abstract

The invention relates to the technical field of biological information, in particular to a methylation sequencing data processing method and system. The method comprises the following steps: acquiring local sequence characteristic data of a DNA sample, wherein the local sequence characteristic data comprises GC content data of 20 basic groups before and after each CpG site, ion concentration data and environment temperature data; performing conformation state prediction on the DNA molecules according to the local sequence characteristic data to obtain thermodynamic parameter matrix data containing melting energy, base accumulation energy and ion acting force; on the basis of the Brown kinetics theory, a motion track model of DNA molecules under thermal disturbance is constructed according to thermodynamic parameter matrix data, and conformation state transition data of all CpG sites are obtained. The sequencing deviation is simulated and corrected through Brown dynamics according to DNA sequence characteristics and thermodynamic parameters, so that the accuracy of methylation sequencing data is improved.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics technology, and in particular to a method and system for processing methylation sequencing data. Background Art

[0002] Methylation sequencing is a technique used to analyze the DNA methylation status. DNA methylation refers to the addition of a methyl (CH3) group at certain positions of DNA, especially on cytosine (C) nucleotides. This process is part of epigenetics and usually affects gene expression without changing the DNA sequence. In methylation sequencing, DNA is usually processed, and certain chemical methods are used to convert unmethylated cytosine into uracil (U), while methylated cytosine does not undergo this conversion. Then, high-throughput sequencing technology is used to sequence the DNA, and by comparing the sequencing results with the reference genome, it is studied which regions of the DNA have undergone methylation changes. The processing process of methylation sequencing data includes multiple steps (such as quality control, methylation data alignment, identification and differential analysis of methylation regions, and calculation of methylation levels, etc.), involving the entire process from the quality control of raw data to biological interpretation.

[0003] However, the traditional methods for processing methylation sequencing data often have the following problems: In whole-genome methylation sequencing, the sequencing depths of different regions may vary greatly. For example, CpG island regions tend to obtain a relatively high coverage, while regions with a low GC content have relatively insufficient coverage. This non-uniformity will lead to inaccurate estimation of methylation levels in some regions. Especially in single-cell sequencing, due to the small amount of starting material, this problem will be more prominent. Summary of the Invention

[0004] Based on this, it is necessary for the present invention to provide a method and system for processing methylation sequencing data to solve at least one of the above technical problems.

[0005] To achieve the above object, a method for processing methylation sequencing data includes the following steps: Step S1: Obtain local sequence feature data of a DNA sample, including GC content data, ion concentration data, and environmental temperature data of 20 bases before and after each CpG site; Step S2: Predict the conformational state of the DNA molecule according to the local sequence feature data to obtain thermodynamic parameter matrix data including melting energy, base stacking energy, and ionic force; Step S3: Based on the Brownian dynamics theory, construct a motion trajectory model of the DNA molecule under thermal perturbation according to the thermodynamic parameter matrix data to obtain conformational state transition data of each CpG site; Step S4: Calculate the theoretical sequencing accessibility index for each CpG site based on the conformational state transition data to obtain sequencing bias prediction data; obtain the actual sequencing depth data, and compare the actual sequencing depth data with the sequencing bias prediction data to obtain depth correction coefficient data; Step S5: Perform adaptive correction based on the methylation level according to the depth correction coefficient data, so as to obtain the methylation level data of CpG sites.

[0006] The present invention provides comprehensive basic information for subsequent analysis by obtaining local sequence feature data of a DNA sample, including the GC content, ion concentration, and environmental temperature data of 20 bases before and after each CpG site. These data can reflect the chemical composition, ion environment, and temperature conditions of the DNA sequence, helping to deeply understand the physicochemical properties of the DNA molecule, and thus providing accurate input parameters for subsequent conformational state prediction and improving the prediction accuracy. Predict the conformational state of the DNA molecule based on the local sequence feature data to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ionic interaction force. This process can quantitatively evaluate the stability of the DNA molecule, reveal its conformational characteristics under different conditions, provide key thermodynamic parameters for subsequent construction of the motion trajectory model, make the model more in line with the actual situation, and enhance the prediction ability of the DNA molecule behavior. Construct a motion trajectory model of the DNA molecule under thermal perturbation based on the Brownian dynamics theory to obtain the conformational state transition data of each CpG site. This model can simulate the dynamic behavior of the DNA molecule under thermal perturbation, reveal the transition law of its conformational state, provide important conformational information for subsequent sequencing bias prediction, and help to more accurately evaluate the possible biases in the sequencing process. Calculate the theoretical sequencing accessibility index for each CpG site based on the conformational state transition data to obtain sequencing bias prediction data, and obtain depth correction coefficient data through comparison with the actual sequencing depth data. This process can combine theoretical prediction with actual sequencing data, accurately identify the biases in the sequencing process, and formulate corresponding correction strategies, improve the accuracy and reliability of the sequencing data, and provide a reliable data basis for subsequent methylation level correction. Perform adaptive correction based on the methylation level according to the depth correction coefficient data to obtain the methylation level data of CpG sites. By comprehensively considering the sequencing bias and methylation level and adopting an adaptive correction method, the methylation state of CpG sites can be more accurately inferred, the accuracy and credibility of the methylation level data can be improved, and a more reliable basis can be provided for subsequent biological research and clinical applications.

[0007] Preferably, step S1 includes the following steps: Step S11: Pretreat the DNA sample and extract the base sequence information to obtain single-stranded DNA sequence data. The pretreatment specifically involves denaturing the double-stranded DNA structure. In the embodiment of the present invention, the double-stranded DNA extracted from the biological sample is placed in a reactor containing Tris-EDTA buffer (pH about 8.0), and an appropriate amount of denaturing agent (such as urea or methanol) is added. Then, it is subjected to high-temperature treatment for 5 minutes in a thermal cycler preheated to 95°C to break the hydrogen bonds between the DNA double strands. Subsequently, the reactor is quickly transferred to an ice bath and cooled for 3 minutes to prevent recombination. Then, the single-stranded DNA is captured and purified using the magnetic bead method (utilizing magnetic particles with a silica coating). Finally, the nucleic acid concentration is measured at a wavelength of 260 nm using a UV spectrophotometer, and the base sequence of the purified single-stranded DNA is read through the Sanger or high-throughput sequencing platform to obtain high-quality single-stranded DNA sequence data.

[0008] Step S12: Calculate the GC content of 20 bases before and after each site in the single-stranded DNA sequence data using the sliding window method to obtain GC content data. In the embodiment of the present invention, after obtaining the single-stranded DNA sequence data, the sliding window algorithm is implemented using a programming language (such as Python combined with the BioPython library). The window size is set to 41 bases (with the center being the target CpG site and 20 bases before and after it). The bases within each window are counted, and the GC content is obtained by calculating the ratio of the number of G and C bases to the total number of bases in the window. If an unknown base N is detected within the window, the program automatically moves the window forward by one base position until there is no N in the window or N is still detected after moving continuously five times. At this time, the GC content value of the last valid window is used to replace the current window data, ultimately ensuring that the GC content data in the neighborhood of each CpG site is both accurate and continuous.

[0009] Step S13: Deploy an ion concentration detector array in the DNA sample buffer, and obtain the spatial distribution of sodium ions, magnesium ions, potassium ions, chloride ions, and calcium ions through real-time monitoring to obtain ion concentration data. In the embodiment of the present invention, a detector array composed of highly sensitive ion-selective electrodes (ISEs) is uniformly deployed in the DNA sample buffer. These electrodes are specifically calibrated for sodium, magnesium, potassium, chloride, calcium and other ions. The potential signals collected in real time are converted into specific concentration values through a pre-established standard curve (describing the linear relationship between the electrode potential and the logarithm of the ion concentration according to the Nernst equation). The sampling frequency is set to 1 Hz, and a microcontroller is used to integrate the data of each electrode in the form of spatial coordinates to generate a detailed data map reflecting the spatial distribution of ions in the sample.

[0010] Step S14: Collect and record the temperature information during DNA sample sequencing to obtain environmental temperature data; In the embodiment of the present invention, a temperature sensor with an accuracy of up to ±0.1 °C (such as a platinum resistance thermometer or a thermocouple) is embedded in the DNA sequencing reactor. The sensor is fixed at multiple key positions in the reactor to comprehensively capture environmental temperature fluctuations. The sensor collects temperature data once per second. After the data is filtered by digital filtering (such as using a moving average filtering algorithm) to remove noise, all temperature readings are with accurate timestamps and are uploaded to the central data processing system in real time, ensuring that the temperature changes during the entire sequencing process are continuously and accurately recorded, providing reliable environmental parameters for subsequent data analysis.

[0011] Step S15: Combine the GC content data, ion concentration data, and environmental temperature data into local sequence feature data.

[0012] In the embodiment of the present invention, the GC content data, ion concentration data, and environmental temperature data obtained respectively in the foregoing steps are subjected to spatio-temporal alignment and normalization processing through data fusion software. The local environment of each CpG site is represented as a feature vector. Among them, the GC content is obtained through the G, C ratio within the window, the ion concentration data is determined according to the average concentration value in the neighborhood of the CpG site corresponding to the position of each electrode, and the temperature data takes the real-time recorded value corresponding to the sequencing time. To ensure the rationality of data fusion, the weighted average method (such as according to the weight coefficients determined by preliminary experiments, the weight of the GC content is set to 0.5, the ion concentration is 0.3, and the temperature is 0.2) is used to normalize and integrate each parameter, and finally a local sequence feature data set formatted as a CSV or JSON file is generated. This data set provides high-precision initial input data for subsequent DNA conformation state prediction.

[0013] The present invention provides a basis for subsequent analysis by preprocessing a DNA sample and extracting base sequence information to obtain single-stranded DNA sequence data. The denaturation treatment dissociates the double-stranded DNA structure, facilitating subsequent extraction of characteristic data from the single-stranded DNA sequence and ensuring the accuracy and reliability of the data. The GC content of 20 bases before and after each site in the single-stranded DNA sequence is calculated using a sliding window method to obtain GC content data. This method can accurately reflect the distribution of the GC content in the DNA sequence, providing important sequence feature information for subsequent conformational state prediction and helping to improve the accuracy of the prediction. An ion concentration detector array is arranged in the DNA sample buffer to monitor the spatial distribution of sodium ions, magnesium ions, potassium ions, chloride ions, and calcium ions in real time, obtaining ion concentration data. This step can accurately obtain the ion concentration information in the environment where the DNA molecule is located, providing key ion force data for subsequent calculation of thermodynamic parameters and helping to more accurately evaluate the stability of the DNA molecule. The temperature information during DNA sample sequencing is collected and recorded to obtain ambient temperature data. Temperature is an important factor affecting the conformation and stability of DNA molecules. Obtaining accurate ambient temperature data can provide a basis for subsequent correction of thermodynamic parameters and ensure the accuracy and reliability of the model. The GC content data, ion concentration data, and ambient temperature data are combined into local sequence feature data, providing comprehensive input information for subsequent conformational state prediction. The integration of these data can more comprehensively reflect the physicochemical properties of DNA molecules, providing a more accurate basis for subsequent analysis and prediction and improving the accuracy and reliability of the entire method.

[0014] Preferably, step S12 includes the following steps: The GC content of 20 bases before and after each site in the single-stranded DNA sequence data is calculated using a sliding window method to obtain GC content data; the specific operation of GC content calculation is that when the 20 bases before and after the site contain an unknown base N, the sliding window is moved forward by 1 position until the sliding window does not contain the unknown base N. If the unknown base N is still contained after moving forward 5 times, the GC content value of the sliding window closest to the upstream of the current site where the GC content calculation has been completed is used as the GC content data of the current site.

[0015] In the embodiments of the present invention, a sliding window method implemented by programming is used to traverse the obtained single-stranded DNA sequence. For each target site, a window with a fixed length of 41 bases (i.e., 20 bases before and after the target site) is constructed centered on it. The total number of G and C bases is counted within this window, and this value is divided by 41 to calculate the GC content of this window, thus forming preliminary GC content data. In specific operations, when one or more unknown bases N are detected within this window, the system will automatically shift the entire window forward by 1 base position, that is, the start and end positions of the window are shifted one base to the right at the same time, and it is rechecked whether the new window still contains the unknown base N. This process will be repeated, and a full-window scan is performed after each shift to ensure no interference from N. If after 5 consecutive shifts, the unknown base N is still detected within the new window, the algorithm will no longer continue to move, but use the GC content data obtained from the previous successfully calculated window without N interference in this region as the GC content of the current site, so as to ensure the continuity and accuracy of the data. In practical applications, this algorithm can be implemented in programming environments such as Python. Using loop structures and string processing functions, the bases within each window are counted one by one, while setting the window size to a fixed 41 bases and the maximum number of consecutive shifts to 5, ensuring that the calculation of the GC content for each CpG site has high robustness and consistency. The obtained GC content data is expressed in percentage or decimal form, providing accurate and reliable basic data for subsequent DNA conformation state prediction and thermodynamic parameter calculation.

[0016] The present invention uses a sliding window method to calculate the GC content of 20 bases before and after each site in the single-stranded DNA sequence data, and can accurately obtain the local GC content information of each site. When encountering an unknown base N, by shifting the sliding window forward by 1 position, the interference of N on the GC content calculation is minimized to ensure the accuracy of the calculation result. Even in the case where the unknown base N is still included after shifting forward 5 times, using the GC content value of the sliding window that has recently completed the GC content calculation upstream of the current site as the GC content data of the current site can also ensure that each site has a reasonable GC content value, avoiding data loss or calculation errors caused by the presence of N. This method can effectively handle the unknown bases in the sequence, ensure the integrity and reliability of the GC content data, provide accurate sequence feature information for subsequent conformation state prediction, and thus improve the accuracy and stability of the entire method.

[0017] The present invention also provides a processing system for methylation sequencing data for performing the above-mentioned processing method for methylation sequencing data. The processing system for methylation sequencing data includes: A sequence feature extraction module, which is used to obtain local sequence feature data of a DNA sample, including GC content data, ion concentration data, and environmental temperature data of 20 bases before and after each CpG site; A conformation prediction module, which is used to predict the conformation state of a DNA molecule according to the local sequence feature data, and obtain thermodynamic parameter matrix data including melting energy, base stacking energy, and ion interaction force; A motion simulation module, which is used to construct a motion trajectory model of a DNA molecule under thermal perturbation based on the Brownian dynamics theory according to the thermodynamic parameter matrix data, and obtain conformation state transition data of each CpG site; A sequencing bias correction module, which is used to calculate the theoretical sequencing accessibility index of each CpG site based on the conformation state transition data to obtain sequencing bias prediction data; obtain actual sequencing depth data, and compare the actual sequencing depth data with the sequencing bias prediction data to obtain depth correction coefficient data; A methylation level correction module, which is used to perform adaptive correction based on the methylation level according to the depth correction coefficient data, so as to obtain CpG site methylation level data.

[0018] The present invention provides comprehensive basic information for subsequent analysis by obtaining local sequence feature data of a DNA sample, including the GC content, ion concentration, and environmental temperature data of 20 bases before and after each CpG site. These data can reflect the chemical composition, ion environment, and temperature conditions of the DNA sequence, facilitating a deeper understanding of the physicochemical properties of DNA molecules. Based on the local sequence feature data, conformational state prediction of the DNA molecule is performed to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ion interaction forces. This process can quantitatively evaluate the stability of the DNA molecule, reveal its conformational characteristics under different conditions, and provide key thermodynamic parameters for the subsequent construction of the motion trajectory model. Based on the Brownian dynamics theory, a motion trajectory model of the DNA molecule under thermal perturbation is constructed according to the thermodynamic parameter matrix data, and conformational state transition data of each CpG site are obtained. This model can simulate the dynamic behavior of the DNA molecule under thermal perturbation, reveal the transition law of its conformational state, and provide important conformational information for the subsequent prediction of sequencing bias. Based on the conformational state transition data, the theoretical sequencing accessibility index of each CpG site is calculated to obtain sequencing bias prediction data, and by comparing with the actual sequencing depth data, depth correction coefficient data are obtained. This process can combine theoretical prediction with actual sequencing data, accurately identify biases in the sequencing process, and formulate corresponding correction strategies to improve the accuracy and reliability of sequencing data. Adaptive correction based on the methylation level is performed according to the depth correction coefficient data to obtain CpG site methylation level data. By comprehensively considering sequencing bias and methylation level and adopting an adaptive correction method, the methylation state of CpG sites can be more accurately inferred, improving the accuracy and credibility of methylation level data and providing a more reliable basis for subsequent biological research and clinical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Other features, objects, and advantages of the present invention will become more apparent by reading the detailed description of non-limiting embodiments with reference to the following drawings: Figure 1 It is a schematic flowchart of the steps of the method for processing methylation sequencing data of the present invention; Figure 2 is Figure 1 a detailed schematic flowchart of step S1 in Figure 3 is Figure 1 a detailed schematic flowchart of step S2 in DETAILED DESCRIPTION OF THE EMBODIMENTS

[0020] The technical method of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention.

[0021] In addition, the accompanying drawings are only schematic diagrams of the present invention and are not necessarily drawn to scale. The same reference numerals in the drawings represent the same or similar parts, and thus repeated descriptions thereof will be omitted. Some of the block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. The functional entities can be implemented in software form, or in one or more hardware modules or integrated circuits, or in different networks and / or processor methods and / or microcontroller methods.

[0022] It should be understood that although terms such as "first" and "second" may be used here to describe various units, these units should not be limited by these terms. These terms are only used to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, the first unit can be called the second unit, and similarly the second unit can be called the first unit. The term "and / or" used here includes any and all combinations of one or more of the listed associated items.

[0023] To achieve the above object, please refer to Figures 1 to 3 , the present invention provides a method for processing methylation sequencing data, and the method includes the following steps: Step S1: Obtain local sequence feature data of a DNA sample, including GC content data, ion concentration data, and ambient temperature data of 20 bases before and after each CpG site. In the embodiment of the present invention, the DNA sample is preprocessed. The double-stranded DNA is denatured by high temperature and chemical reagents to obtain a single-stranded DNA sequence; then a high-throughput sequencer is used to obtain accurate base sequence information, and the sliding window method (the window size is fixed at 41 bases, that is, 20 bases on each side of the target CpG site) is used to calculate the GC content in each window. When there is an unknown base N in the window, the system will automatically shift the window forward by one position until it is shifted forward five times continuously and there is still N, then the GC content value of the previous window is used instead. At the same time, an ion concentration detector array is arranged in the DNA sample buffer to monitor and record the concentrations of sodium, magnesium, potassium, chlorine, calcium and other ions in the sample in real time, and a high-precision temperature sensor is built into the sample reactor to obtain ambient temperature data. Finally, the above GC content, ion concentration and temperature data are integrated to form local sequence feature data for subsequent processing.

[0024] Step S2: Predict the conformational state of the DNA molecule based on the local sequence feature data to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ionic interaction force; In the embodiment of the present invention, based on the local sequence feature data obtained in step S1, first, by reading the GC content data of the region where each CpG site is located, the base pair pairing energy of this region is estimated using a method based on the thermodynamic free energy model. Subsequently, in combination with the ion concentration data, the electrolyte distribution model is used to correct the DNA double-strand stability. Then, the nearest-neighbor thermodynamic parameter method is used to calculate the stacking interaction between adjacent bases, thereby obtaining the base stacking energy data. In addition, the ionic interaction force is adjusted for temperature dependence based on the environmental temperature data. The integrated melting energy, base stacking energy, and ionic interaction force information form a multi-dimensional thermodynamic parameter matrix, which can comprehensively reflect the conformational energy distribution of the DNA molecule under experimental conditions.

[0025] Step S3: Based on the Brownian dynamics theory, construct a motion trajectory model of the DNA molecule under thermal perturbation according to the thermodynamic parameter matrix data to obtain the conformational state transition data of each CpG site; In the embodiment of the present invention, using the thermodynamic parameter matrix constructed in step S2, first, the DNA molecular backbone is discretized into a flexible chain model composed of multiple nodes and the elastic connections between them. Then, according to the Brownian dynamics theory, thermal noise is introduced. The random motion of each node under thermal perturbation is described by constructing the motion equation of the node. The velocity Verlet algorithm is used for time integration to update the node positions. At the same time, the interaction forces generated by bond length stretching, bond angle bending, and torsional effects between adjacent nodes are calculated within each time step to simulate the motion trajectory of the entire DNA molecule under specific ion concentration and temperature conditions. Furthermore, by statistically analyzing the conformational changes corresponding to the CpG sites within each time step, the conformational state transition data is obtained.

[0026] Step S4: Calculate the theoretical sequencing accessibility index of each CpG site based on the conformational state transition data to obtain sequencing bias prediction data; obtain the actual sequencing depth data, and compare the actual sequencing depth data with the sequencing bias prediction data to obtain depth correction coefficient data; Based on the conformational state transition data obtained in step S3, the embodiments of the present invention first quantitatively evaluate the local exposure degree and unwinding characteristics of the DNA double strand in different conformational states, construct a theoretical sequencing accessibility index, that is, map the empirical correspondence between the local conformational characteristics of DNA and the sequencing reaction efficiency, and obtain the prediction of the sequencing accessibility of each CpG site under ideal conditions; then, obtain the actual sequencing depth data through high-throughput sequencing technology, calculate the ratio of the actual depth to the predicted theoretical index, set the correction coefficient to 1 when the ratio is between 0.8 and 1.2, and use the reciprocal of the ratio as the correction coefficient when the ratio exceeds this range, and finally obtain a depth correction coefficient data that can reflect the sequencing deviation.

[0027] Step S5: Perform adaptive correction based on the methylation level according to the depth correction coefficient data, so as to obtain the methylation level data of CpG sites.

[0028] The embodiments of the present invention use the depth correction coefficient data obtained in step S4 to preliminarily count the methylation signals and non-methylation signals in the original methylation sequencing reads, and use Bayesian statistical methods to comprehensively consider the sequencing depth and signal ratio to obtain the preliminary methylation level data of CpG sites; subsequently, in combination with the context environment and conformational characteristics of the local DNA sequence, identify and correct the outliers in the preliminary data, group the CpG sites with similar conformational characteristics by cluster analysis, calculate the regional correction parameters, and then perform adaptive refinement correction on the methylation level of each CpG site according to these regional parameters, and finally form accurate methylation level data that can truly reflect the methylation state of the DNA sample.

[0029] The present invention provides comprehensive basic information for subsequent analysis by obtaining local sequence feature data of a DNA sample, including the GC content, ion concentration, and environmental temperature data of 20 bases before and after each CpG site. These data can reflect the chemical composition, ion environment, and temperature conditions of the DNA sequence, helping to deeply understand the physicochemical properties of the DNA molecule, and thus providing accurate input parameters for subsequent conformational state prediction and improving the prediction accuracy. Based on the local sequence feature data, the conformational state of the DNA molecule is predicted to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ion interaction force. This process can quantitatively evaluate the stability of the DNA molecule, reveal its conformational characteristics under different conditions, provide key thermodynamic parameters for the subsequent construction of the movement trajectory model, make the model more in line with the actual situation, and enhance the prediction ability of the DNA molecule behavior. Based on the Brownian dynamics theory, a movement trajectory model of the DNA molecule under thermal perturbation is constructed to obtain the conformational state transition data of each CpG site. This model can simulate the dynamic behavior of the DNA molecule under thermal perturbation, reveal the transition law of its conformational state, provide important conformational information for subsequent sequencing bias prediction, and help to more accurately evaluate the possible biases in the sequencing process. Based on the conformational state transition data, the theoretical sequencing accessibility index of each CpG site is calculated to obtain sequencing bias prediction data, and through comparison with the actual sequencing depth data, depth correction coefficient data are obtained. This process can combine theoretical prediction with actual sequencing data, accurately identify the biases in the sequencing process, and formulate corresponding correction strategies to improve the accuracy and reliability of the sequencing data, providing a reliable data basis for subsequent methylation level correction. According to the depth correction coefficient data, an adaptive correction based on the methylation level is performed to obtain the methylation level data of the CpG site. By comprehensively considering the sequencing bias and methylation level and adopting an adaptive correction method, the methylation state of the CpG site can be more accurately inferred, improving the accuracy and credibility of the methylation level data, and providing a more reliable basis for subsequent biological research and clinical applications.

[0030] Preferably, step S1 includes the following steps: Step S11: Pretreat the DNA sample and extract the base sequence information to obtain single-stranded DNA sequence data, where the pretreatment is specifically to denature the double-stranded DNA structure; In an embodiment of the present invention, double-stranded DNA extracted from a biological sample is placed in a reactor containing Tris-EDTA buffer (pH about 8.0), and an appropriate amount of denaturing agent (such as urea or methanol) is added. Then, it is subjected to high-temperature treatment in a thermal cycler preheated to 95°C for 5 minutes to break the hydrogen bonds between the DNA double strands. Subsequently, the reactor is quickly transferred to an ice bath and cooled for 3 minutes to prevent recombination. Next, the magnetic bead method (using magnetic particles with a silica coating) is used to capture and purify single-stranded DNA. Finally, a UV spectrophotometer is used to measure the nucleic acid concentration at a wavelength of 260 nm, and the base sequence of the purified single-stranded DNA is read through a Sanger or high-throughput sequencing platform, thereby obtaining high-quality single-stranded DNA sequence data.

[0031] Step S12: Calculate the GC content of 20 bases before and after each site in the single-stranded DNA sequence data using the sliding window method, thereby obtaining GC content data. In an embodiment of the present invention, after obtaining the single-stranded DNA sequence data, the sliding window algorithm is implemented using a programming language (such as Python combined with the BioPython library). The window size is set to 41 bases (with the center being the target CpG site and 20 bases before and after it). The bases within each window are counted, and the GC content is obtained by calculating the ratio of the number of G and C bases to the total number of bases in the window. If an unknown base N is detected within the window, the program automatically moves the window forward by one base position until there is no N in the window or N is still detected after moving continuously five times. At this time, the GC content value of the last valid window is used to replace the current window data, ultimately ensuring that the GC content data in the neighborhood of each CpG site is both accurate and continuous.

[0032] Step S13: Arrange an ion concentration detector array in the DNA sample buffer, and obtain the spatial distribution of sodium ions, magnesium ions, potassium ions, chloride ions, and calcium ions through real-time monitoring, thereby obtaining ion concentration data. In an embodiment of the present invention, a detector array composed of high-sensitivity ion-selective electrodes (ISEs) is uniformly arranged in the DNA sample buffer. These electrodes are specifically calibrated for sodium, magnesium, potassium, chloride, calcium and other ions. Through a pre-established standard curve (describing the linear relationship between the electrode potential and the logarithm of the ion concentration according to the Nernst equation), the collected potential signals are converted into specific concentration values in real time. The sampling frequency is set to 1 Hz, and a microcontroller is used to integrate the data of each electrode in the form of spatial coordinates to generate a detailed data map reflecting the spatial distribution of ions in the sample.

[0033] Step S14: Collect and record the temperature information during DNA sample sequencing, thereby obtaining ambient temperature data. In the embodiments of the present invention, a temperature sensor with an accuracy of up to ±0.1 °C (such as a platinum resistance thermometer or a thermocouple) is embedded in the DNA sequencing reactor. The sensor is fixed at multiple key positions in the reactor to comprehensively capture environmental temperature fluctuations. The sensor collects temperature data once per second. After the data is filtered by digital filtering (for example, using a moving average filtering algorithm) to remove noise, all temperature readings are with accurate timestamps and are uploaded to the central data processing system in real time, ensuring that the temperature changes during the entire sequencing process are continuously and accurately recorded, providing reliable environmental parameters for subsequent data analysis.

[0034] Step S15: Combine the GC content data, ion concentration data, and environmental temperature data into local sequence feature data.

[0035] In the embodiments of the present invention, the GC content data, ion concentration data, and environmental temperature data respectively obtained in the foregoing steps are subjected to spatio-temporal alignment and normalization processing by data fusion software. The local environment of each CpG site is represented as a feature vector, where the GC content is obtained through the G, C ratio within the window, the ion concentration data is determined according to the average concentration value in the neighborhood of the CpG site corresponding to the position of each electrode, and the temperature data takes the real-time recorded value corresponding to the sequencing time. To ensure the rationality of data fusion, the weighted average method (for example, according to the weight coefficients determined by preliminary experiments, the weight of the GC content is set to 0.5, the ion concentration is 0.3, and the temperature is 0.2) is used to normalize and integrate each parameter, and finally a local sequence feature data set formatted as a CSV or JSON file is generated. This data set provides high-precision initial input data for subsequent DNA conformation state prediction.

[0036] The present invention provides a basis for subsequent analysis by preprocessing a DNA sample and extracting base sequence information to obtain single-stranded DNA sequence data. The denaturation treatment dissociates the double-stranded DNA structure, facilitating subsequent extraction of characteristic data from the single-stranded DNA sequence and ensuring the accuracy and reliability of the data. The GC content of 20 bases before and after each site in the single-stranded DNA sequence is calculated using a sliding window method to obtain GC content data. This method can accurately reflect the distribution of the GC content in the DNA sequence, providing important sequence characteristic information for subsequent conformational state prediction and helping to improve the accuracy of the prediction. An ion concentration detector array is arranged in the DNA sample buffer to monitor the spatial distribution of sodium ions, magnesium ions, potassium ions, chloride ions, and calcium ions in real time, obtaining ion concentration data. This step can accurately obtain the ion concentration information in the environment where the DNA molecule is located, providing key ion force data for subsequent calculation of thermodynamic parameters and helping to more accurately evaluate the stability of the DNA molecule. The temperature information during DNA sample sequencing is collected and recorded to obtain ambient temperature data. Temperature is an important factor affecting the conformation and stability of DNA molecules. Obtaining accurate ambient temperature data can provide a basis for subsequent correction of thermodynamic parameters and ensure the accuracy and reliability of the model. The GC content data, ion concentration data, and ambient temperature data are combined into local sequence characteristic data, providing comprehensive input information for subsequent conformational state prediction. The integration of these data can more comprehensively reflect the physicochemical properties of DNA molecules, providing a more accurate basis for subsequent analysis and prediction and improving the accuracy and reliability of the entire method.

[0037] Preferably, step S12 includes the following steps: The GC content of 20 bases before and after each site in the single-stranded DNA sequence data is calculated using a sliding window method to obtain GC content data; specifically, when the 20 bases before and after the site contain an unknown base N, the sliding window is moved forward by 1 position until the sliding window does not contain the unknown base N. If the unknown base N is still contained after moving forward 5 times, the GC content value of the sliding window closest to the upstream of the current site where the GC content calculation has been completed is used as the GC content data of the current site.

[0038] In the embodiments of the present invention, a sliding window method implemented by programming is used to traverse the obtained single-stranded DNA sequence. For each target site, a window with a fixed length of 41 bases (i.e., 20 bases before and after the target site) is constructed centered on it. The total number of G and C bases is counted within this window, and this value is divided by 41 to calculate the GC content of this window, thereby forming preliminary GC content data. In specific operations, when one or more unknown bases N are detected within the window, the system will automatically shift the entire window forward by 1 base position, that is, the start and end positions of the window are both shifted one base to the right, and the new window is rechecked to see if it still contains the unknown base N. This process will be repeated, and a full-window scan is performed after each shift to ensure that there is no interference from N. If after 5 consecutive shifts, the unknown base N is still detected in the new window, the algorithm will no longer continue to move, but use the GC content data obtained from the previous successfully calculated window without N interference in this region as the GC content of the current site, so as to ensure the continuity and accuracy of the data. In practical applications, this algorithm can be implemented in programming environments such as Python. The bases within each window are counted one by one using loop structures and string processing functions, while setting the window size to a fixed 41 bases and the maximum number of consecutive shifts to 5, ensuring that the GC content calculation for each CpG site has high robustness and consistency. The obtained GC content data is expressed in percentage or decimal form, providing accurate and reliable basic data for subsequent DNA conformation state prediction and thermodynamic parameter calculation.

[0039] The present invention uses a sliding window method to calculate the GC content of 20 bases before and after each site in the single-stranded DNA sequence data, and can accurately obtain the local GC content information of each site. When encountering an unknown base N, by shifting the sliding window forward by 1 position, the interference of N on the GC content calculation is avoided as much as possible to ensure the accuracy of the calculation result. Even in the case where the unknown base N is still included after shifting forward 5 times, using the GC content value of the sliding window that has completed the GC content calculation closest to the upstream of the current site as the GC content data of the current site can also ensure that each site has a reasonable GC content value, avoiding data loss or calculation errors caused by the existence of N. This method can effectively handle the unknown bases in the sequence, ensure the integrity and reliability of the GC content data, provide accurate sequence feature information for subsequent conformation state prediction, and thus improve the accuracy and stability of the entire method.

[0040] Preferably, step S2 includes the following steps: Step S21: Calculate the base pair pairing energy of the region where each CpG site is located according to the local sequence feature data, and quantitatively evaluate the stability of double-stranded DNA through a thermodynamic free energy model, thereby obtaining base pairing energy data; In the embodiments of the present invention, the base sequence information around each CpG site is extracted from the local sequence feature data. Using the existing thermodynamic free energy model, the hydrogen bond stability of each base pair is converted into a quantitative description. That is, by comparing the differences between adjacent nucleotide pairs (such as A-T and G-C), the free energy values measured in the laboratory (for example, the G-C pair usually has a lower free energy, indicating higher stability) are used to evaluate the stability of the double-stranded DNA in the target region. During the whole process, the system gradually scans the base pair combinations in the target region, accumulates the energy contributed by each base pair according to the preset parameters (such as the standard free energy value at 25 °C), and finally outputs a set of base pairing energy data reflecting the overall stability of this region. This data not only considers the stability of a single base pair but also incorporates the modulation effect of the local sequence context on the pairing energy.

[0041] Step S22: Based on the base pairing energy data, perform an analysis of the stacking interaction between adjacent bases, and use the nearest neighbor thermodynamic parameter method to calculate the interaction strength between bases, so as to obtain the base stacking energy data; After obtaining the base pairing energy data in step S21 in the embodiments of the present invention, the system uses the nearest neighbor thermodynamic parameter method to conduct a detailed analysis of the stacking interaction between adjacent bases. This method is based on the experimentally measured data and believes that the - interaction between two consecutive bases makes a significant contribution to the overall stability of DNA. The software will calculate the stacking energy contribution for each pair of adjacent bases according to the pre-established parameter library. This process involves quantifying the stacking effect of each pair of bases through the "adjacent effect" described in words. For example, a higher stacking energy value may be given to a specific sequence combination, while a lower value may be given to other combinations. Finally, the stacking energies of all adjacent bases in the target region are summed up to generate a set of detailed base stacking energy data, reflecting the subtle differences in conformational stability of the local sequence.

[0042] Step S23: Establish an electrolyte distribution model according to the ion concentration data, and perform an evaluation of the ion atmosphere distribution around the DNA molecule based on the Poisson-Boltzmann equation, so as to obtain the electrostatic potential energy data; In an embodiment of the present invention, based on the ion concentration data collected in real time, an electrolyte distribution model is first established. This model uses the ion diffusion coefficients determined in advance in the laboratory and the ion concentrations in the buffer solution, and places the DNA molecule in a hypothetical three-dimensional space grid. It is explained in a text description that: at each grid point, according to the concentrations of sodium, magnesium, potassium, chlorine, and calcium ions measured nearby, the system simulates an ion distribution field, and then uses a numerical solution method based on the Poisson-Boltzmann equation to evaluate the potential distribution around the DNA molecule. Considering the characteristic that ions decrease with distance in the electric field, the system continuously adjusts the potential value of each grid point through an iterative algorithm until the entire scenario reaches an equilibrium state, thereby outputting a set of detailed electrostatic potential energy data, which reflects the influence of the ion distribution around the DNA molecule on the molecular stability.

[0043] Step S24: Perform temperature-dependent correction on the base pairing energy data, base stacking energy data, and electrostatic potential energy data based on the environmental temperature data, so as to obtain temperature-corrected energy parameter data; In an embodiment of the present invention, the environmental temperature data is combined with the previously obtained base pairing energy, base stacking energy, and electrostatic potential energy data. First, the system detects the temperature recorded in real time during the sequencing process and compares this temperature with the laboratory standard temperature, and then performs energy correction according to the temperature-dependent law described in the text. For example, it is explained that an increase in temperature will cause the hydrogen bond vibration to intensify, thereby weakening the pairing energy, or the temperature change will affect the ion shielding effect. The software multiplies each energy data item by a temperature correction factor (these factors are determined by experimental data, such as through Arrhenius-type adjustment or description based on the Debye-Hückel theory), so that the final energy parameter data can accurately reflect the energy state of the DNA molecule under the actual temperature conditions, and the corrected data is more in line with the molecular behavior in a specific experimental environment.

[0044] Step S25: Map the energy parameter data into a three-dimensional space coordinate system to establish a conformational energy landscape of the DNA molecule, so as to obtain conformational energy distribution data; In an embodiment of the present invention, after obtaining the temperature-corrected energy parameter data, the system maps these energy data into a three-dimensional coordinate system through computer graphics tools. This process presents the physical structure of the DNA molecule in the form of spatial coordinates, and each coordinate point corresponds to the energy value at a specific position on the DNA molecule. The system uses an interpolation method to smooth the discrete energy data into a continuous energy surface, forming a conformational energy landscape similar to a topographic map, where the low-energy region represents a more stable conformational state and the high-energy region represents an unstable state. Through this three-dimensional mapping, not only can the energy distribution of the DNA molecule under specific conditions be intuitively displayed, but also it can provide a visual basis and detailed energy distribution data for subsequent conformational prediction.

[0045] Step S26: Sample and predict the possible conformational states of the DNA molecule based on the conformational energy distribution data, so as to obtain the thermodynamic parameter matrix data.

[0046] Based on the conformational energy landscape constructed in step S25, the embodiment of the present invention systematically explores the possible conformational states of the DNA molecule by using a sampling algorithm. The specific method adopts the Metropolis sampling strategy in Monte Carlo simulation, that is, in each iteration, a possible conformational change is simulated by locally perturbing the DNA backbone, and the energy difference before and after this change is calculated. If the energy of the new conformation is lower, it is accepted with a high probability; if the energy increases, it is accepted with a lower probability according to the Boltzmann distribution. During the whole process, the system continuously updates the dihedral angle and bond angle parameters of the DNA, records the energy state of each conformational transition. After a large number of samplings, the system can construct a detailed thermodynamic parameter matrix data. This matrix not only describes the possible conformational states near each CpG site, but also reflects the transition probabilities between these states, thus providing a theoretical basis and numerical support for understanding the dynamic behavior of DNA under thermal perturbation.

[0047] The present invention calculates the base pair pairing energy of the region where each CpG site is located, and uses the thermodynamic free energy model to quantitatively evaluate the stability of double-stranded DNA, obtaining base pairing energy data, which provides key energy parameters for subsequent analysis and helps to accurately evaluate the stability of DNA molecules. Based on the base pairing energy data, the stacking interaction strength between adjacent bases is calculated by the nearest neighbor thermodynamic parameter method, obtaining base stacking energy data, which further enriches the energy parameters of DNA molecules and helps to more comprehensively understand the physicochemical properties of DNA molecules. An electrolyte distribution model is established according to the ion concentration data, and the Poisson-Boltzmann equation is used to evaluate the ion atmosphere distribution around the DNA molecule, obtaining electrostatic potential energy data, which provides information on ion forces for the energy parameters of DNA molecules and helps to more accurately evaluate the behavior of DNA molecules in different ion environments. Based on the environmental temperature data, the base pairing energy data, base stacking energy data, and electrostatic potential energy data are corrected for temperature dependence, obtaining temperature-corrected energy parameter data, ensuring the accuracy and reliability of the energy parameters under different temperature conditions and providing precise data support for the subsequent establishment of the conformational energy landscape. The energy parameter data is mapped into a three-dimensional space coordinate system to establish the conformational energy landscape of the DNA molecule, obtaining conformational energy distribution data, which provides an intuitive visualization basis for predicting the conformational states of DNA molecules and helps to more accurately identify and predict the possible conformational states of DNA molecules. According to the conformational energy distribution data, the possible conformational states of the DNA molecule are sampled and predicted, obtaining thermodynamic parameter matrix data, which provides key thermodynamic parameters for the subsequent construction of the motion trajectory model and helps to more accurately simulate the behavior of DNA molecules under thermal perturbation, thereby improving the accuracy and reliability of the entire method.

[0048] Preferably, step S24 includes the following steps: Step S241: Calculate the degree of influence of temperature change on molecular motion according to the environmental temperature data through the Arrhenius equation, thereby obtaining temperature influence factor data; In the embodiment of the present invention, the real-time environmental temperature data in the reactor is collected, converted into Kelvin temperature, and compared with a preset standard temperature (such as 298 Kelvin). The Arrhenius equation is used to describe the influence of temperature on the molecular motion rate. Through the predetermined activation energy and frequency factor, the system uses exponential decay to describe the trend of the molecular motion rate accelerating when the temperature increases, thereby calculating a temperature influence factor, which reflects the change amplitude of the molecular motion rate relative to the standard state at the current environmental temperature and is represented by a dimensionless value, providing a quantitative basis for the temperature correction of subsequent energy parameters.

[0049] Step S242: Analyze the change in hydrogen bond strength at different temperatures for the base pairing energy data based on the Gibbs free energy equation, so as to obtain the corrected base pairing energy data; In the embodiment of the present invention, starting from the base pairing energy data obtained in step S21, combined with the real-time environmental temperature data, according to the principle of the influence of temperature on the free energy of reaction in the Gibbs free energy equation, the system first conducts a temperature-dependent analysis of the hydrogen bond stability of each base pair. Considering the phenomenon that the probability of hydrogen bond breakage increases and the stability decreases when the temperature rises, by comparing the equilibrium state of hydrogen bond breakage and formation at the current temperature and the standard temperature, the system adjusts the original base pairing energy according to the temperature deviation and outputs the base pairing energy data corrected by temperature. This data more accurately reflects the stability of the DNA double strand at the actual experimental temperature.

[0050] Step S243: Dynamically adjust the base stacking energy data based on the temperature dependence of the stacking interaction between aromatic rings according to the temperature influence factor data, so as to obtain the corrected base stacking energy data; - In the embodiment of the present invention, using the temperature influence factor data obtained in step S241, this step conducts a temperature-dependent adjustment for the π-π stacking interaction between aromatic rings. The system first analyzes the contribution of the aromatic ring stacking interaction to the molecular stability at the standard temperature, and then multiplies the original base stacking energy data by a dynamic adjustment factor according to the molecular motion intensification effect described by the temperature influence factor. This factor will reduce the stacking energy value when the temperature is higher than the standard and increase the stacking energy value when the temperature is lower than the standard, so as to generate the corrected base stacking energy data reflecting the temperature-dependent change, ensuring that this data can accurately reflect the actual change of the interaction force between aromatic rings in a high-temperature or low-temperature environment.

[0051] Step S244: Calculate the influence of temperature change on the ionic atmosphere distribution based on the electrostatic potential energy data according to the Debye-Hückel theory, so as to obtain the corrected electrostatic potential energy data; In the embodiment of the present invention, based on the electrostatic potential energy data obtained in step S23, in accordance with the Debye-Hückel theory, by collecting the temperature data to calculate the thermal motion effect of ions in the solution, the system first determines the ionic shielding length at the standard temperature, and then adjusts the dielectric constant and the ionic activity coefficient according to the principle that the enhanced ionic thermal motion and the weakened shielding effect caused by the increase in temperature, so as to recalculate the local potential distribution and obtain the electrostatic potential energy data corrected by temperature. This data reflects the actual influence of the ionic atmosphere on the electrostatic interaction of DNA under the current temperature conditions.

[0052] ​Step S245: Perform weighted integration processing on the base pairing energy correction data, base stacking energy correction data, and electrostatic potential energy correction data to obtain energy weight data, where the weighted integration processing specifically uses the entropy weight method to determine the weight coefficients of each energy term; In the embodiment of the present invention, after obtaining the base pairing energy, base stacking energy, and electrostatic potential energy data corrected by temperature, the system uses the entropy weight method to perform weighted integration processing on each energy term. First, standardize each group of data, calculate the probability distribution of each data item in the entire sample and obtain its information entropy. According to the principle that data items with lower information entropy have higher information content, assign higher weights to them, while assign lower weights to those with higher information entropy. Finally, integrate each energy term through weighted average to generate a set of energy weight data, which comprehensively reflects the contributions of various energy effects to the DNA stability under temperature changes.

[0053] Step S246: Normalize the base pairing energy correction data, base stacking energy correction data, and electrostatic potential energy correction data according to the energy weight data, and establish a unified temperature correction standard to obtain temperature-corrected energy parameter data.

[0054] In the embodiment of the present invention, the energy weight data obtained in step S245 is used to normalize each corrected energy data. The specific operations include converting each group of data to the standard interval of 0 to 1 by subtracting the minimum value and then dividing by the range, and then performing weighted synthesis according to their respective weights to establish a unified temperature correction standard, thereby outputting a set of standardized temperature-corrected energy parameter data, which provides unified and accurate input parameters for subsequent DNA molecular dynamics simulations and conformation predictions.

[0055] The present invention calculates the influence degree of temperature change on molecular motion according to the Arrhenius equation based on environmental temperature data, obtains temperature influence factor data, provides a basis for subsequent energy parameter correction, and ensures the accuracy and reliability of energy parameters under different temperature conditions. The change analysis of hydrogen bond strength at different temperatures for base pairing energy data is carried out based on the Gibbs free energy equation, and base pairing energy correction data are obtained, further improving the accuracy of base pairing energy data and making it more in line with the actual temperature conditions. The base stacking energy data are dynamically adjusted based on the temperature dependence of π-π stacking between aromatic rings according to the temperature influence factor data, and base stacking energy correction data are obtained, enabling the base stacking energy data to more accurately reflect the influence of temperature on the stacking effect. Based on the Debye-Hückel theory, the influence of temperature change on the ionic atmosphere distribution is calculated according to the electrostatic potential energy data, and electrostatic potential energy correction data are obtained, ensuring the accuracy of electrostatic potential energy data under different temperature conditions and providing reliable data support for subsequent energy parameter integration. The base pairing energy correction data, base stacking energy correction data, and electrostatic potential energy correction data are weighted and integrated to obtain energy weight data. The weighted integration process specifically uses the entropy weight method to determine the weight coefficients of each energy term, ensuring a more reasonable weight distribution of different energy terms during the integration process and improving the comprehensive accuracy of energy parameters. According to the energy weight data, the base pairing energy correction data, base stacking energy correction data, and electrostatic potential energy correction data are normalized to establish a unified temperature correction standard, and temperature-corrected energy parameter data are obtained, ensuring the comparability and consistency of different energy parameters after temperature correction, providing accurate energy parameters for subsequent construction of conformational energy landscapes, and thus improving the accuracy and reliability of the entire method.

[0056] Preferably, step S26 includes the following steps: Step S261: Perform an energy gradient analysis on the conformational energy distribution data, and determine the set of stable conformational states of the DNA molecule based on the principle of physical and chemical energy minimization, so as to obtain the initial conformational state data; In the embodiment of the present invention, a detailed energy gradient analysis is performed on the discrete energy data obtained from the conformational energy landscape. The numerical differentiation method is used to calculate the change rate of energy with spatial position in each local region. The system quickly locates the energy minimum point through the gradient descent method, and identifies local energy flat regions within the region where the energy gradient tends to zero. These regions are regarded as the set of possible stable conformational states; in actual operation, set the threshold of the energy gradient (for example, less than kcal / ) to ensure that the selected conformation is indeed in the energy minimum state, and thus output the initial conformational state data, representing the most stable conformation of the DNA molecule under specific experimental conditions.

[0057] Step S262: Design a Metropolis sampling strategy based on the initial conformation state data, and perform small deformations on the DNA molecular backbone through a local perturbation algorithm to obtain conformation perturbation sequence data; Based on the initial conformation state data obtained in Step S261, the embodiment of the present invention systematically designs a Metropolis sampling strategy, locally perturbs the DNA molecular backbone within each sampling period. The perturbation operations mainly involve fine-tuning local dihedral angles, bond angles, and a small number of bond lengths, and the perturbation amplitude is strictly controlled within the range of ±3 to ±5 degrees to ensure the smallness and continuity of the structural deformation. After each perturbation, the differences between the new conformation state and the original state are recorded to form a series of continuous conformation perturbation sequence data, providing preliminary data support for subsequent energy change calculations and conformation conversions.

[0058] Step S263: Calculate the energy change of conformation conversion based on the conformation perturbation sequence data, and use the Boltzmann distribution to determine the acceptance probability of the new conformation, thereby obtaining conformation conversion probability data; In the embodiment of the present invention, using the conformation perturbation sequence data obtained in Step S262, the system calculates the energy difference between the front and rear conformations after each conformation change. By substituting the energy change into the Boltzmann distribution formula (where the transition probability is proportional to exp(- ), k is the Boltzmann constant, and T is the absolute temperature), the acceptance probability of the new conformation after each perturbation is determined, thereby establishing a conformation conversion probability data set, which reflects the thermodynamic feasibility and stability of each step of conformation transition.

[0059] Step S264: Perform Markov chain iterative sampling according to the conformation conversion probability data, and update the dihedral angle and bond angle parameters of the DNA molecule in each sampling step, thereby obtaining conformation evolution trajectory data; After obtaining the conformation conversion probability data in the embodiment of the present invention, the system constructs a Markov chain model, and updates the key conformation parameters of the DNA molecule (mainly including dihedral angles and bond angles) through iterative sampling. Each step of the update is based on the conversion probability calculated in the previous step. If the new conformation is accepted with a high probability, the current state is updated; otherwise, the original state is maintained. After a large number of iterative samplings, the system records the trajectory of each step of conformation change to form detailed conformation evolution trajectory data, comprehensively reflecting the dynamic conformation evolution process of the DNA molecule under thermal perturbation conditions.

[0060] Step S265: Perform the main conformation conversion path based on the principal component analysis method on the conformation evolution trajectory data to obtain conformation conversion path data; In the embodiment of the present invention, for the high-dimensional conformational evolution trajectory data obtained in step S264, the system uses the principal component analysis (PCA) method for dimensionality reduction. By calculating the variance contribution of each variable in the data, several principal components with the highest percentage of total variation are extracted, and then the most important conversion paths in the conformational change are identified. In specific operations, a cumulative contribution rate threshold (such as reaching more than 80%) is set to ensure that the extracted main paths can accurately represent the overall conformational transition trend, and finally one or more key conformational conversion path data are obtained, intuitively showing the transition process of the DNA molecule between different energy states.

[0061] Step S266: Based on the conformational conversion path data, construct a thermodynamic parameter matrix including the energy difference, conversion probability, and conformational stability index between adjacent conformational states, so as to obtain the thermodynamic parameter matrix data.

[0062] In the embodiment of the present invention, based on the conformational conversion path data obtained in step S265, the system integrates the energy difference between adjacent conformational states, the conversion probability obtained through Boltzmann distribution and Markov chain sampling, and the stability index of each conformational state (such as local energy minimum, energy gradient flatness, etc.) to construct a multi-dimensional thermodynamic parameter matrix. In the specific implementation process, each parameter is first normalized, and then the weighted average method (the weights are determined in advance according to the importance of each parameter) is used to fuse each data item into the matrix, so that each matrix element not only reflects the energy barrier between adjacent conformational states, but also reflects the conversion possibility and stability index, thereby generating a thermodynamic parameter matrix data that comprehensively describes the conformational conversion behavior of the DNA molecule under thermal perturbation, providing a detailed theoretical basis for subsequent conformational prediction and kinetic simulation.

[0063] The present invention performs energy gradient analysis on conformational energy distribution data, and determines the set of stable conformational states of a DNA molecule based on the principle of physicochemical energy minimization to obtain initial conformational state data, providing an accurate starting point for subsequent conformational sampling and prediction, and ensuring the rationality and effectiveness of the sampling process. According to the initial conformational state data, a Metropolis sampling strategy is designed, and the DNA molecular backbone is slightly deformed through a local perturbation algorithm to obtain conformational perturbation sequence data, which can effectively explore its conformational space while maintaining the rationality of the DNA molecular structure, increasing the diversity and comprehensiveness of the sampling process. The energy change during conformational conversion is calculated based on the conformational perturbation sequence data, and the acceptance probability of the new conformation is determined using the Boltzmann distribution to obtain conformational conversion probability data, making the conformational conversion process conform to the principles of thermodynamics and ensuring the physical authenticity and reliability of the sampling process. Markov chain iterative sampling is performed according to the conformational conversion probability data, and the dihedral angle and bond angle parameters of the DNA molecule are updated in each sampling step to obtain conformational evolution trajectory data, which can record in detail the evolution process of the DNA molecular conformation and provide rich data support for subsequent conformational analysis. The main conformational conversion path analysis based on the principal component analysis method is performed on the conformational evolution trajectory data to obtain conformational conversion path data, which can effectively extract the main features and trends during the conformational conversion process, simplify complex data, and highlight key information. Based on the conformational conversion path data, a thermodynamic parameter matrix containing the energy difference, conversion probability, and conformational stability index between adjacent conformational states is constructed to obtain thermodynamic parameter matrix data, providing comprehensive and accurate thermodynamic parameters for the conformational state prediction of the DNA molecule, making the prediction results more accurate and reliable, and providing a solid foundation for the subsequent construction of the motion trajectory model.

[0064] Preferably, step S3 includes the following steps: Step S31: According to the thermodynamic parameter matrix data, the DNA backbone is discretized into a flexible chain composed of nodes and elastic connections to construct DNA backbone model data; In the embodiment of the present invention, first, the continuous backbone of the DNA molecule is discretized into a series of finite nodes according to the thermodynamic parameter matrix data. Each node represents a base pair or a nucleotide sequence fragment of a fixed length. Then, elastic connections are introduced between adjacent nodes. The equilibrium lengths and stiffness coefficients of these virtual springs are preset according to the experimentally determined energy parameters (for example, the distance between each node can be set to about 3.4 angstroms, and the stiffness coefficient is adjusted according to the free energy data), thereby constructing discrete backbone model data that reflects the local conformation and overall flexibility characteristics of the DNA molecule, providing a basic structure for subsequent dynamic simulations.

[0065] Step S32: Based on the DNA backbone model data, establish the Brownian dynamics equation of motion according to the motion characteristics of each point under thermodynamic fluctuations, so as to obtain the node motion equation data; In the embodiment of the present invention, the DNA backbone model data obtained in S31 is utilized. The system establishes the motion equation of each node according to the Brownian dynamics theory. Assuming that the node is simultaneously driven by the elastic restoring force and the random thermal noise, the Langevin equation containing the friction force and the random force term is introduced to describe the motion characteristics of the node, and the corresponding parameters are set according to the temperature and medium viscosity in the experimental environment (for example, the friction coefficient and the random force amplitude are set). Finally, the motion equation data describing the motion behavior of each node is generated, providing a theoretical basis for simulating the random motion of DNA under thermal perturbation.

[0066] Step S33: Calculate the interaction forces between the nodes according to the node motion equation data, and describe them with the many-body potential energy function, so as to obtain the node force data, where the interaction forces include bond length stretching force, bond angle bending force and torsional force; In the embodiment of the present invention, based on the node motion equation data obtained in S32, the system further calculates the interaction forces between adjacent nodes, mainly considering three aspects: First, the bond length stretching force is evaluated by measuring the deviation between the distance between nodes and the preset equilibrium distance; Second, the bond angle bending force is determined by calculating the deviation between the angle formed by three consecutive nodes and the ideal angle; Finally, the torsional force is estimated by the change in the dihedral angle determined by four consecutive nodes. All these forces are described by the many-body potential energy function, and its core idea is to quantify the proportional relationship between each force term and the square of the deviation. After numerical integration calculation, a set of complete node force data is formed, reflecting the comprehensive effect of the interaction between each node in the DNA backbone.

[0067] Step S34: Use the velocity Verlet algorithm to perform time evolution on the motion equation according to the node force data, and simulate the motion trajectory of the DNA molecule under thermal perturbation, so as to obtain the molecular motion trajectory data; In the embodiment of the present invention, after obtaining the node force data, the system uses the velocity Verlet algorithm to perform time stepping solution on the motion equation of each node. This algorithm can update the position and velocity of the node simultaneously and maintain a high energy conservation accuracy. In the specific operation, the time step is set to about 0.1 femtosecond. By calculating the acceleration of the current node according to the force condition at each time step, and updating the position and velocity in turn. After a large number of iterations, the system simulates the continuous motion trajectory of the DNA molecule under thermal perturbation conditions, and finally outputs a series of detailed molecular motion trajectory data, recording the spatial position changes of each node at different time points.

[0068] Step S35: Extract the conformational features from the molecular motion trajectory data. By calculating the time evolution of local conformational parameters, the conformational feature sequence data is obtained, where the local conformational parameters include the bending angle, the torsion angle, and the persistence length. Based on the molecular motion trajectory data obtained in S34, the embodiments of the present invention adopt a data processing algorithm to extract the local conformational parameters of the DNA backbone at each time step. The local conformational parameters include the bending angle calculated from three consecutive nodes, the torsion angle determined by four consecutive nodes, and the persistence length describing the length of the local straight segment. The sliding window and time series statistical methods are used to gradually calculate the evolution trend of these parameters over time. Finally, these local conformational information is integrated into a set of conformational feature sequence data, providing a quantitative description for revealing the regularity of DNA molecular conformational changes.

[0069] Step S36: Identify the main conformational states of the DNA molecule and their transition characteristics according to the conformational feature sequence data. Analyze the kinetic process of conformational transition by using the Markov state model, so as to obtain the conformational state transition data.

[0070] The embodiments of the present invention utilize the conformational feature sequence data extracted in S35. The system first identifies a series of representative stable conformational states through clustering analysis and data segmentation methods, and then establishes a Markov state model. Assuming that the DNA conformational transition process conforms to the Markov property, by statistically analyzing the transition frequencies between stable states in consecutive time steps, a state transition matrix is constructed, where each matrix element reflects the probability of transitioning from one stable state to another. After optimizing the model parameters using the maximum likelihood estimation, a set of conformational state transition data is finally output. This data details the main conformational states of the DNA molecule under thermal perturbation and their transition dynamics, providing a solid theoretical support for subsequent conformational prediction and kinetic simulation.

[0071] Based on the thermodynamic parameter matrix data, the present invention constructs the DNA backbone model data by discretizing the DNA backbone into a flexible chain composed of nodes and elastic connections, providing a basic framework for subsequent motion simulation and being able to accurately describe the physical structure and elastic properties of DNA molecules. Based on the DNA backbone model data, the Brownian dynamics motion equation is established to obtain the node motion equation data, which can accurately describe the motion characteristics of DNA molecules under the action of thermodynamic fluctuations and provide a theoretical basis for simulating the motion trajectory of DNA molecules. Calculate the interaction forces between each node and describe them with a many-body potential energy function to obtain the node force data, including bond length stretching force, bond angle bending force, and torsional force, which can comprehensively consider various interaction forces inside DNA molecules and ensure the accuracy and reliability of motion simulation. Using the velocity Verlet algorithm, the motion equation is time-evolved according to the node force data to simulate the motion trajectory of DNA molecules under thermal perturbation and obtain the molecular motion trajectory data, which can efficiently perform long-term motion simulation and capture the dynamic behavior of DNA molecules under thermal perturbation. Extract the conformational characteristic data from the molecular motion trajectory data. By calculating the time evolution of local conformational parameters, the conformational characteristic sequence data, including bending angle, torsional angle, and persistence length, can be obtained, which can describe in detail the local conformational changes of DNA molecules and provide key information for subsequent conformational state recognition. According to the conformational characteristic sequence data, identify the main conformational states of DNA molecules and their conversion characteristics, and analyze the kinetic process of conformational conversion through a Markov state model to obtain the conformational state transition data, which can accurately reveal the conversion rules between different conformational states of DNA molecules and provide important conformational information for subsequent sequencing bias prediction.

[0072] Preferably, step S4 includes the following steps: Step S41: Analyze the exposure degree and local unwinding characteristics of the DNA double strand in different conformational states according to the conformational state transition data, so as to obtain the conformational accessibility data; In the embodiment of the present invention, the time series information extracted from the conformational state transition data is used to analyze the spatial configuration of the local region of the DNA double strand in each conformational state, and a pre-established judgment criterion is adopted to evaluate the exposure degree of the double strand, that is, by calculating the theoretical solvent accessible area and local unwinding angle of each region (for example, when the local torsional angle exceeds a certain threshold, such as 20 degrees, it is regarded as a partially unwound state), and combining the opening and closing states of base pairs recorded in the molecular dynamics simulation, the exposure and opening degrees of key regions (especially near CpG sites) in each conformational state are statistically analyzed. Finally, these quantitative indicators are normalized and output as conformational accessibility data, providing a direct structural basis for subsequent sequencing efficiency prediction.

[0073] Step S42: Calculate the theoretical sequencing accessibility index for each CpG site based on the conformational accessibility data, thereby obtaining the theoretical accessibility index data, where the theoretical sequencing accessibility index specifically extracts the correspondence between the conformational state and the sequencing efficiency; In the embodiment of the present invention, based on the conformational accessibility data obtained in step S41, the system constructs a mapping model to correspond the accessibility of the region where each CpG site is located to the expected sequencing efficiency. The specific operation is to utilize the statistical relationship between the sequencing efficiency and the local conformational state in the known experimental data (such as using linear regression or non - linear fitting methods) to determine the theoretical sequencing accessibility index for each CpG site. This index is usually represented by a value between 0 and 1, where a higher value indicates that the site is more likely to be captured by the sequencer in theory. Finally, a set of theoretical accessibility index data covering all CpG sites is obtained.

[0074] Step S43: Predict the sequencing bias by combining the local conformational characteristics of the DNA molecule with the sequencing reaction kinetics according to the theoretical accessibility index data, thereby obtaining the sequencing bias prediction data; In the embodiment of the present invention, after obtaining the theoretical accessibility index data, the system comprehensively analyzes the local conformational characteristics of each CpG site (such as exposure degree, unwinding state, and local geometric parameters) and the sequencing reaction kinetics parameters (such as amplification rate, reaction time, and enzyme activity). Through multivariate regression or machine learning model training, the possible deviation situations of each CpG site in the actual sequencing reaction are predicted. The specific method is to use the theoretical accessibility index as the main input variable and combine other kinetic parameters for fitting. Finally, the sequencing bias prediction data corresponding to each CpG site is output, and this data reflects the influence trend of the conformational state on the sequencing efficiency.

[0075] Step S44: Perform high - throughput sequencing on the pre - obtained sequencing library, and perform filtering and pre - processing to obtain the actual sequencing depth data; In the embodiment of the present invention, a high - throughput sequencing platform (such as the Illumina NovaSeq or HiSeq system) is used to perform parallel sequencing on the pre - constructed sequencing library. After the raw data is generated, it is filtered and pre - processed through standard bioinformatics processes, including removing low - quality sequences, removing adapters and duplicate reads, and performing genome alignment and correction processing. Finally, the actual sequencing coverage data for each CpG site is statistically obtained. This data is presented in the form of read counts for each site and stored in a CSV or other standard data format to ensure the accuracy and reproducibility of the data.

[0076] Step S45: Calculate the ratio of the actual sequencing depth data to the sequencing bias prediction data to obtain the sequencing bias ratio data; In the embodiment of the present invention, the actual sequencing depth data obtained in step S44 is used to calculate the ratio with the sequencing deviation data predicted in step S43. The specific method is to calculate the ratio of the actual sequencing depth to the theoretical predicted value for each CpG site respectively. This ratio reflects the deviation degree between the conformational influence and the sequencing reaction in actual operation. If the ratio is close to 1, it indicates that the sequencing depth is consistent with the theoretical prediction. If the ratio is significantly greater than or less than 1, it indicates that there is a deviation, providing a quantitative basis for further correction, and finally generating a set of complete sequencing deviation ratio data.

[0077] Step S46: Formulate an adaptive correction strategy according to the sequencing deviation ratio data, set different correction coefficients for different deviation ranges, so as to obtain the depth correction coefficient data. The adaptive correction strategy is specifically as follows: when the sequencing deviation ratio is greater than 1.2, the correction coefficient is the reciprocal of the sequencing deviation ratio; when the sequencing deviation ratio is less than 0.8, the correction coefficient is the reciprocal of the sequencing deviation ratio; when the sequencing deviation ratio is between 0.8 and 1.2, the correction coefficient is 1.

[0078] Based on the sequencing deviation ratio data obtained in step S45, the embodiment of the present invention systematically formulates an adaptive correction strategy, that is, classifies the deviation ratios of each CpG site. When the ratio is greater than 1.2, it is considered that the actual sequencing depth is higher than the theoretical expectation. At this time, the correction coefficient is set to the reciprocal of this ratio to reduce the depth value; when the ratio is less than 0.8, it indicates that the actual depth is lower than the expectation, and the reciprocal strategy is also adopted to increase the depth value; when the ratio is between 0.8 and 1.2, the correction coefficient is set to 1, indicating that no adjustment is required. The specific parameters (such as 1.2 and 0.8 are empirical thresholds) are determined through the training and verification of a large number of standard sample data. Finally, a set of depth correction coefficient data for each CpG site is output, providing an accurate quantitative tool for the subsequent fine correction of the methylation level.

[0079] The present invention analyzes the exposure degree and local unwinding characteristics of DNA double strands in different conformational states based on conformational state transition data, obtains conformational accessibility data, provides key conformational information for subsequent sequencing bias prediction, and can accurately evaluate the accessibility of DNA molecules in different conformational states. Based on the conformational accessibility data, the theoretical sequencing accessibility index of each CpG site is calculated to obtain theoretical accessibility index data. By extracting the corresponding relationship between the conformational state and sequencing efficiency, a theoretical basis for sequencing bias prediction is provided, making the prediction results more accurate and reliable. According to the theoretical accessibility index data, the local conformational characteristics of DNA molecules are combined with the kinetics of the sequencing reaction to predict sequencing bias, obtaining sequencing bias prediction data, which can comprehensively consider the conformational characteristics of DNA molecules and the kinetic factors of the sequencing reaction, and improve the accuracy and reliability of sequencing bias prediction. The pre-obtained sequencing library is subjected to high-throughput sequencing, filtered and preprocessed to obtain actual sequencing depth data, providing actual sequencing data for subsequent calculation of the sequencing bias ratio, and ensuring the accuracy and reliability of the data. The actual sequencing depth data and the sequencing bias prediction data are used for ratio calculation to obtain sequencing bias ratio data, which can intuitively reflect the degree of bias existing in the sequencing process and provide a basis for formulating subsequent adaptive correction strategies. An adaptive correction strategy is formulated according to the sequencing bias ratio data, and different correction coefficients are set for different bias ranges to obtain depth correction coefficient data. Through specific correction strategies, sequencing bias can be effectively corrected, the accuracy and reliability of sequencing data can be improved, and a solid foundation is provided for subsequent methylation level correction.

[0080] Preferably, step S5 includes the following steps: Step S51: Perform linear correction processing on the sequencing bias prediction data according to the depth correction coefficient data to obtain sequencing depth data; In the embodiment of the present invention, the sequencing bias ratio data obtained from step S45 and the previously calculated depth correction coefficient data are acquired, and a linear correction method is used for each CpG site, that is, the theoretical sequencing bias prediction value of the site is multiplied by the corresponding correction coefficient to obtain the corrected sequencing depth data; for example, when the bias prediction value of a certain site is 120 reads and the correction coefficient is 0.9, the system calculates their product as 108 reads as the actual sequencing depth of the site. This process is performed in the software in a batch matrix operation manner to ensure rapid and accurate linear adjustment of data for tens of thousands or even hundreds of thousands of CpG site data, and all data are stored in a standardized format for subsequent analysis.

[0081] Step S52: Obtain the methylation and non-methylation signal counts in the original methylation sequencing reads to obtain the original methylation signal data; In the embodiment of the present invention, by performing data parsing and filtering on the original methylation sequencing reads, the reads are screened using preset quality control parameters (such as setting the minimum sequencing quality value to 20 and the mapping quality not less than 30). Then, the position of each CpG site is located through a comparison software, and the base information of each site appearing in the original data is counted. The count of maintaining C (representing the methylation signal) and the count of being converted to T (representing the non-methylation signal) at the corresponding position are respectively statistically analyzed. Finally, these count results are organized into a structured original data table of methylation signals, with each row corresponding to a CpG site, facilitating subsequent statistical analysis and Bayesian inference processing.

[0082] Step S53: Based on the sequencing depth data and the original methylation signal data, perform methylation status inference based on Bayesian statistics to obtain preliminary methylation level data. In the embodiment of the present invention, using the sequencing depth data obtained in step S51 and the original methylation and non-methylation signal data obtained in step S52, the system adopts the Bayesian statistical method for methylation status inference. In this process, the observed data of each CpG site is regarded as success and failure events in a binomial distribution experiment. A suitable beta distribution is selected as the prior distribution, and the posterior distribution is obtained by updating the prior information. Then, the mean value of the posterior distribution is used as the preliminary methylation level data of this site. The whole process fully considers sequencing errors and sample variability, and the data update and probability calculation are automatically realized through an iterative algorithm in the software to ensure that the methylation probability value of each CpG site accurately reflects its true biological state.

[0083] Step S54: Perform local sequence context correlation analysis on the preliminary methylation level data, identify and correct methylation level outliers affected by DNA conformation to obtain denoised methylation data. In the embodiment of the present invention, based on the preliminary methylation level data obtained in step S53, correlation analysis is performed in combination with the context information of the local sequence where each CpG site is located (such as the surrounding GC content, DNA conformation state, etc.). Outliers that deviate significantly from the average level of adjacent regions are identified through a sliding window and local regression analysis method. A deviation threshold (for example, more than or less than 0.15 from the local average) is set as the abnormal determination criterion. Once abnormal data is detected, the system uses the median of adjacent CpG sites or re-estimates and corrects through a local smoothing algorithm. Finally, a denoised methylation data after noise removal and abnormal correction is output, laying a foundation for further regional correction.

[0084] Step S55: Based on the denoised methylation data, perform CpG site clustering analysis to identify CpG groups with similar conformation characteristics, thereby obtaining regional correction parameter data. After obtaining the denoised methylation data in the embodiments of the present invention, the system uses a clustering analysis method to group all CpG sites. Specifically, in the operation, the K-means or hierarchical clustering algorithm is used to perform multi-dimensional feature clustering according to the methylation level of each site and its corresponding DNA conformation features (such as the degree of local unwinding, exposure index, etc.), and a predetermined number of clusters is set or the optimal number of clusters is automatically determined using the silhouette coefficient, so as to identify CpG groups with similar conformation features and methylation patterns. The features within each group are highly consistent. Subsequently, regional correction parameters (such as mean, standard deviation, and correction coefficient) are calculated for each group to perform overall correction on the data within the group and improve the signal confidence in subsequent processing.

[0085] Step S56: Perform fine correction of the signal confidence on the denoised methylation data according to the regional correction parameter data, so as to obtain the methylation level data of CpG sites.

[0086] In the embodiments of the present invention, according to the regional correction parameter data determined in step S55, fine signal confidence correction processing is performed on the denoised methylation data of each CpG site. In this process, a weighted adjustment method is used to combine the original methylation level of each CpG site with the correction parameters of its region, and the value is re-corrected in a manner such as weighted average or linear regression model. The consistency of the methylation level within the region and the stability of the local conformation features are considered during the correction process. The output result is the final methylation level data of CpG sites. These data not only reflect the methylation status of individual sites but also take into account the overall biological characteristics of the region, providing more accurate and reliable epigenetic information.

[0087] The present invention performs linear correction processing on the sequencing bias prediction data according to the depth correction coefficient data to obtain the sequencing depth data, which can effectively correct the sequencing bias, improve the accuracy and reliability of the sequencing data, and provide a solid foundation for subsequent methylation level inference. The methylation and non-methylation signal counts in the original methylation sequencing reads are obtained to get the original methylation signal data, which provides the original data support for subsequent methylation state inference and ensures the integrity and reliability of the data. According to the sequencing depth data and the original methylation signal data, methylation state inference based on Bayesian statistics is performed to obtain the preliminary methylation level data. By comprehensively considering the sequencing depth and methylation signal, the accuracy and reliability of methylation state inference are improved using the Bayesian statistical method. Local sequence context correlation analysis is performed on the preliminary methylation level data to identify and correct the methylation level outliers affected by DNA conformation, and the denoised methylation data is obtained, which can effectively remove the methylation level outliers caused by DNA conformation and improve the accuracy and credibility of the methylation level data. CpG site clustering analysis is performed based on the denoised methylation data to identify CpG groups with similar conformational characteristics, and the regional correction parameter data is obtained. By clustering analysis, CpG groups with similar conformational characteristics are identified, which provides parameter support for subsequent regional correction and improves the pertinence and effectiveness of the correction. According to the regional correction parameter data, fine correction of the signal confidence of the denoised methylation data is performed to obtain the CpG site methylation level data. By fine correction, the accuracy and reliability of the methylation level data are further improved, providing high-quality data support for the final methylation level analysis.

[0088] The present invention also provides a processing system for methylation sequencing data for performing the above-mentioned processing method of methylation sequencing data. The processing system for methylation sequencing data includes: A sequence feature extraction module for obtaining local sequence feature data of a DNA sample, including GC content data, ion concentration data, and environmental temperature data of 20 bases before and after each CpG site; A conformation prediction module for predicting the conformation state of a DNA molecule according to the local sequence feature data to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ion interaction force; A motion simulation module for constructing a motion trajectory model of a DNA molecule under thermal perturbation based on the Brownian dynamics theory according to the thermodynamic parameter matrix data to obtain the conformation state transition data of each CpG site; A sequencing bias correction module for calculating the theoretical sequencing accessibility index of each CpG site based on the conformation state transition data to obtain sequencing bias prediction data; obtaining the actual sequencing depth data, and comparing the actual sequencing depth data with the sequencing bias prediction data to obtain depth correction coefficient data; The methylation level correction module is used to perform adaptive correction based on the methylation level according to the depth correction coefficient data, so as to obtain the methylation level data of CpG sites.

[0089] In the present invention, by obtaining the local sequence feature data of a DNA sample, including the GC content, ion concentration, and ambient temperature data of 20 bases before and after each CpG site, comprehensive basic information is provided for subsequent analysis. These data can reflect the chemical composition, ion environment, and temperature conditions of the DNA sequence, and contribute to an in-depth understanding of the physicochemical properties of DNA molecules. Based on the local sequence feature data, the conformational state of the DNA molecule is predicted to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ion interaction force. This process can quantitatively evaluate the stability of the DNA molecule, reveal its conformational characteristics under different conditions, and provide key thermodynamic parameters for the subsequent construction of the motion trajectory model. Based on the Brownian dynamics theory, according to the thermodynamic parameter matrix data, a motion trajectory model of the DNA molecule under thermal perturbation is constructed to obtain the conformational state transition data of each CpG site. This model can simulate the dynamic behavior of the DNA molecule under thermal perturbation, reveal the transition law of its conformational state, and provide important conformational information for the subsequent prediction of sequencing bias. Based on the conformational state transition data, the theoretical sequencing accessibility index of each CpG site is calculated to obtain the sequencing bias prediction data, and through comparison with the actual sequencing depth data, the depth correction coefficient data is obtained. This process can combine theoretical prediction with actual sequencing data, accurately identify the bias in the sequencing process, and formulate corresponding correction strategies to improve the accuracy and reliability of sequencing data. According to the depth correction coefficient data, adaptive correction based on the methylation level is performed to obtain the methylation level data of CpG sites. By comprehensively considering the sequencing bias and methylation level and adopting an adaptive correction method, the methylation state of CpG sites can be inferred more accurately, improving the accuracy and credibility of the methylation level data, and providing a more reliable basis for subsequent biological research and clinical applications.

[0090] Therefore, from any perspective, the embodiments should be regarded as exemplary and non-limiting. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, all changes falling within the meaning and scope of the equivalent elements of the application documents are intended to be encompassed within the present invention.

[0091] The above description is only the specific implementation manners of the present invention, enabling those skilled in the art to understand or implement the present invention. Various modifications to these embodiments will be obvious to those skilled in the art. The general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features invented herein.

Claims

1. A method for processing methylation sequencing data, characterized in that, Including the following steps: Step S1: Obtain the local sequence feature data of the DNA sample, including the GC content data, ion concentration data, and environmental temperature data of 20 bases before and after each CpG site; Step S2: Predict the conformational state of the DNA molecule based on the local sequence feature data to obtain the thermodynamic parameter matrix data including the melting energy, base stacking energy, and ionic interaction force; Step S3: Based on the Brownian dynamics theory, construct the motion trajectory model of the DNA molecule under thermal perturbation according to the thermodynamic parameter matrix data to obtain the conformational state transition data of each CpG site; Step S4: Calculate the theoretical sequencing accessibility index of each CpG site based on the conformational state transition data to obtain the sequencing bias prediction data; Obtain the actual sequencing depth data, and compare the actual sequencing depth data with the sequencing bias prediction data to obtain the depth correction coefficient data; Step S5: Perform adaptive correction based on the methylation level according to the depth correction coefficient data to obtain the methylation level data of the CpG site.

2. The method for processing methylation sequencing data according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Pretreat the DNA sample and extract the base sequence information to obtain the single-stranded DNA sequence data, where the pretreatment is specifically to denature the double-stranded DNA structure; Step S12: Use the sliding window method to calculate the GC content of 20 bases before and after each site in the single-stranded DNA sequence data to obtain the GC content data; Step S13: Arrange an ion concentration detector array in the DNA sample buffer, and obtain the spatial distribution of sodium ions, magnesium ions, potassium ions, chloride ions, and calcium ions through real-time monitoring to obtain the ion concentration data; Step S14: Collect and record the temperature information during DNA sample sequencing to obtain the environmental temperature data; Step S15: Combine the GC content data, ion concentration data, and environmental temperature data into local sequence feature data.

3. The method for processing methylation sequencing data according to claim 2, characterized in that, Step S12 includes the following steps: Use the sliding window method to calculate the GC content of 20 bases before and after each site in the single-stranded DNA sequence data to obtain the GC content data; The specific operation of the GC content calculation is that when there is an unknown base N among the 20 bases before and after the site, move the sliding window forward by 1 position until there is no unknown base N in the sliding window. If there is still an unknown base N after moving forward 5 times, then use the GC content value of the sliding window closest to the upstream of the current site where the GC content calculation has been completed as the GC content data of the current site.

4. The method for processing methylation sequencing data according to claim 3, characterized in that, Step S2 includes the following steps: Step S21: Calculate the base pair pairing energy of the region where each CpG site is located according to the local sequence feature data, and quantitatively evaluate the stability of double-stranded DNA through the thermodynamic free energy model to obtain the base pairing energy data; Step S22: Perform the stacking interaction analysis between adjacent bases based on the base pairing energy data, and calculate the interaction strength between bases using the nearest neighbor thermodynamic parameter method to obtain the base stacking energy data; Step S23: Establish an electrolyte distribution model based on the ion concentration data, and evaluate the ion atmosphere distribution around the DNA molecule based on the Poisson-Boltzmann equation, so as to obtain the electrostatic potential energy data; Step S24: Perform temperature-dependent correction on the base pairing energy data, base stacking energy data, and electrostatic potential energy data based on the environmental temperature data, so as to obtain the temperature-corrected energy parameter data; Step S25: Map the energy parameter data into a three-dimensional space coordinate system, establish the conformational energy landscape of the DNA molecule, so as to obtain the conformational energy distribution data; Step S26: Sample and predict the possible conformational states of the DNA molecule according to the conformational energy distribution data, so as to obtain the thermodynamic parameter matrix data.

5. The method for processing methylation sequencing data according to claim 4, characterized in that, Step S24 includes the following steps: Step S241: Calculate the influence degree of temperature change on molecular motion according to the environmental temperature data through the Arrhenius equation, so as to obtain the temperature influence factor data; Step S242: Analyze the change of hydrogen bond strength at different temperatures for the base pairing energy data based on the Gibbs free energy equation, so as to obtain the base pairing energy correction data; Step S243: Based on the temperature influence factor data, the base stacking energy data is subjected to the temperature influence factor data. - Dynamic adjustment of the temperature dependence of the stacking effect, resulting in base stacking energy correction data; Step S244: Calculate the influence of temperature change on the ion atmosphere distribution based on the Debye-Hückel theory according to the electrostatic potential energy data, so as to obtain the electrostatic potential energy correction data; Step S245: Perform weighted integration processing on the base pairing energy correction data, base stacking energy correction data, and electrostatic potential energy correction data, so as to obtain the energy weight data, where the weighted integration processing is specifically to determine the weight coefficients of each energy term by using the entropy weight method; Step S246: Perform normalization processing on the base pairing energy correction data, base stacking energy correction data, and electrostatic potential energy correction data according to the energy weight data, establish a unified temperature correction standard, so as to obtain the temperature-corrected energy parameter data.

6. The method for processing methylation sequencing data according to claim 5, characterized in that, Step S26 includes the following steps: Step S261: Perform energy gradient analysis on the conformational energy distribution data, and determine the set of stable conformational states of the DNA molecule based on the principle of minimizing physical and chemical energy, so as to obtain the initial conformational state data; Step S262: Design a Metropolis sampling strategy according to the initial conformational state data, and perform small deformations on the DNA molecule backbone through a local perturbation algorithm, so as to obtain the conformational perturbation sequence data; Step S263: Calculate the energy change of conformational conversion based on the conformational perturbation sequence data, and determine the acceptance probability of the new conformation by using the Boltzmann distribution, so as to obtain the conformational conversion probability data; Step S264: Perform Markov chain iterative sampling according to the conformational conversion probability data, and update the dihedral angle and bond angle parameters of the DNA molecule in each sampling step, so as to obtain the conformational evolution trajectory data; Step S265: Perform the main conformational conversion path based on the principal component analysis method on the conformational evolution trajectory data, so as to obtain the conformational conversion path data; Step S266: Construct a thermodynamic parameter matrix containing the energy differences, transition probabilities, and conformational stability indices between adjacent conformational states based on the conformational transition path data, thereby obtaining the thermodynamic parameter matrix data.

7. The processing method of methylation sequencing data according to claim 6, wherein, Step S3 includes the following steps: Step S31: Based on the thermodynamic parameter matrix data, construct the DNA backbone model data by discretizing the DNA backbone into a flexible chain composed of nodes and elastic connections. Step S32: Based on the DNA backbone model data, establish the Brownian dynamics equation of motion according to the motion characteristics of each point under thermodynamic fluctuations, thereby obtaining the node motion equation data. Step S33: Calculate the interaction forces between the nodes according to the node motion equation data and describe them with a many-body potential function, thereby obtaining the node force data, where the interaction forces include bond length stretching force, bond angle bending force, and torsional force. Step S34: Use the velocity Verlet algorithm to perform time evolution on the equation of motion according to the node force data, and simulate the motion trajectory of the DNA molecule under thermal perturbation, thereby obtaining the molecular motion trajectory data. Step S35: Extract the conformational characteristics from the molecular motion trajectory data, and obtain the conformational characteristic sequence data by calculating the time evolution of the local conformational parameters, where the local conformational parameters include bending angle, torsional angle, and persistence length. Step S36: Identify the main conformational states and their transition characteristics of the DNA molecule according to the conformational characteristic sequence data, and analyze the kinetic process of conformational transition by the Markov state model, thereby obtaining the conformational state transition data.

8. The processing method of methylation sequencing data according to claim 7, wherein, Step S4 includes the following steps: Step S41: Analyze the exposure degree and local unwinding characteristics of the DNA double strand in different conformational states according to the conformational state transition data, thereby obtaining the conformational accessibility data. Step S42: Calculate the theoretical sequencing accessibility index of each CpG site based on the conformational accessibility data, thereby obtaining the theoretical accessibility index data, where the theoretical sequencing accessibility index specifically extracts the corresponding relationship between the conformational state and the sequencing efficiency. Step S43: Predict the sequencing bias by combining the local conformational characteristics of the DNA molecule with the sequencing reaction kinetics according to the theoretical accessibility index data, thereby obtaining the sequencing bias prediction data. Step S44: Perform high-throughput sequencing on the pre-acquired sequencing library, and perform filtering and preprocessing, thereby obtaining the actual sequencing depth data. Step S45: Calculate the ratio of the actual sequencing depth data to the sequencing bias prediction data, thereby obtaining the sequencing bias ratio data. Step S46: Develop an adaptive correction strategy according to the sequencing bias ratio data, and set different correction coefficients for different bias ranges, thereby obtaining the depth correction coefficient data, where the adaptive correction strategy is specifically: when the sequencing bias ratio is greater than 1.2, the correction coefficient is the reciprocal of the sequencing bias ratio; when the sequencing bias ratio is less than 0.8, the correction coefficient is the reciprocal of the sequencing bias ratio; when the sequencing bias ratio is between 0.8 and 1.2, the correction coefficient is 1.

9. The processing method of methylation sequencing data according to claim 8, wherein, Step S5 includes the following steps: Step S51: Perform linear correction processing on the sequencing bias prediction data according to the depth correction coefficient data to obtain the sequencing depth data; Step S52: Obtain the methylation and non-methylation signal counts in the original methylation sequencing reads to obtain the original methylation signal data; Step S53: Infer the methylation status based on Bayesian statistics according to the sequencing depth data and the original methylation signal data to obtain the preliminary methylation level data; Step S54: Perform local sequence context correlation analysis on the preliminary methylation level data, identify and correct the methylation level outliers affected by DNA conformation to obtain the denoised methylation data; Step S55: Perform CpG site clustering analysis based on the denoised methylation data, identify CpG groups with similar conformational characteristics to obtain the regional correction parameter data; Step S56: Perform fine correction of the signal confidence on the denoised methylation data according to the regional correction parameter data to obtain the CpG site methylation level data.

10. A processing system for methylation sequencing data, wherein, A processing system for performing the processing method of methylation sequencing data as described in claim 1, the processing system for methylation sequencing data includes: A sequence feature extraction module, configured to obtain local sequence feature data of a DNA sample, including GC content data, ion concentration data, and environmental temperature data of 20 bases before and after each CpG site; A conformation prediction module, configured to predict the conformation state of a DNA molecule according to the local sequence feature data to obtain thermodynamic parameter matrix data including the melting energy, base stacking energy, and ionic force; A motion simulation module, configured to construct a motion trajectory model of a DNA molecule under thermal perturbation based on the Brownian dynamics theory according to the thermodynamic parameter matrix data to obtain the conformation state transition data of each CpG site; A sequencing bias correction module, configured to calculate the theoretical sequencing accessibility index of each CpG site based on the conformation state transition data to obtain the sequencing bias prediction data; obtain the actual sequencing depth data, and compare the actual sequencing depth data with the sequencing bias prediction data to obtain the depth correction coefficient data; A methylation level correction module, configured to perform adaptive correction based on the methylation level according to the depth correction coefficient data to obtain the CpG site methylation level data.

Images