Variation recognition method for nanopore single molecule sequencing current signals

Abstract

This invention relates to a method for identifying variants in current signals from nanopore single-molecule sequencing. The method divides the raw current signal generated by nanopore single-molecule sequencing into continuous time periods according to the order in which bases pass through the nanopore. The average signal rhythm is calculated, and a standard time rhythm is constructed for the standard base sequence passing through the nanopore. The deviation difference between the measured average signal rhythm and the standard time rhythm is compared; regions with significant differences are identified as potential variant regions. Current signals from the preceding and following neighborhoods are extracted, and the perturbation difference is quantified. Regions with significant directional asymmetry in the perturbation difference are identified as variant regions to be confirmed. The perturbation intensity is calculated by expanding the current signal extraction range. When the perturbation intensity shows a continuous trend of gradually increasing or decreasing intensity, the variant region to be confirmed is deemed physically plausible. The theoretical current signal is determined using the corresponding unmutated sequence. The measured current signal is compared, and the current difference is calculated. Regions exceeding the normal difference range are confirmed as true variant regions.

Description

Technical Field

[0001] This invention relates to nanopore single-molecule sequencing current signals, specifically a method for identifying variations in nanopore single-molecule sequencing current signals. Background Technology

[0002] Existing technologies, such as the nanopore sequencing signal evaluation method, device, electronic device and storage medium proposed in Chinese patent CN117594130A, have certain signal quality evaluation and process monitoring capabilities, but their technical paths and structures have obvious shortcomings and limitations in adaptability, making it difficult to meet the requirements for highly sensitive, robust and spatially accurate identification of current signal variations. First, this invention focuses on segmenting and evaluating nanopore sequencing signals from a macroscopic process perspective. Its method divides the raw current signal into blank current segments, bound current segments, and continuous through-pore current segments, and performs statistical analysis on the signal fluctuations, mean, and rate of change of these different states to extract evaluation indicators such as average current value, standard deviation, and settling time. This allows the method to determine whether the signal at a certain time meets the sequencing requirements or whether there is an equipment malfunction. This type of method essentially serves sequencing quality assessment and instrument status judgment, rather than being designed for fine-grained variation detection at the base level. Therefore, it lacks structural support in identifying variations such as small insertions, deletions, and single base substitutions. Its indicator dimensions cannot cover current behaviors that are highly correlated with mutations, such as the changing patterns of current rhythm and the spatial propagation characteristics of structural perturbations.

[0003] Secondly, the statistical indicators used in this method are based on coarse-grained statistical characteristics of the entire signal, focusing particularly on state identification before and after nanopore passage, rather than dynamically tracking base channel perturbations during continuous passage. While this approach is suitable for macroscopic state judgment under instrument calibration or low-pass filtering conditions, it neglects the fine-grained perturbation characteristics of the current signal and fails to form a systematic analytical framework for the intensity of perturbations in the current and their performance at different window scales. For example, in high-offset perturbation signals, which are typically manifested as single or a few consecutive sampling points, it is difficult to capture local deviations and perturbation trend changes using only the average and variance statistics of the original signal. Therefore, it is impossible to analyze and judge signal indicators reflecting changes in sequence structure, such as rhythm shifts and rate anomalies, resulting in insufficient recognition ability or even false alarms when facing complex regions such as high GC regions, repetitive sequences, and structurally variable regions.

[0004] Further analysis reveals that the technical solution fails to incorporate multi-scale analysis mechanisms, making it impossible to uniformly determine the distribution patterns of anomalous structures at different granularities. In current sequencing processes, base mutations and structural anomalies often manifest simultaneously as local perturbations and cross-regional trend shifts. This solution lacks cross-scale correlation analysis methods and does not establish correlation modeling for perturbation intensity spectra at different window scales. Consequently, when variant sites are close together, morphologically inconsistent, or have significant cumulative effects, boundary blurring or misclassification can easily occur, reducing actual localization accuracy. Furthermore, its method for identifying perturbation signal boundaries only relies on statistically significant numerical interval division, lacking support based on perturbation intensity trend changes, sequence continuity analysis, and change point detection algorithms. This prevents precise determination of the start and end points of variant regions, thus limiting its application in base localization, variant length estimation, and reliability ranking. Finally, from an engineering practicality perspective, while emphasizing sequencing signal quality control, this solution does not construct a complete data flow loop around variant identification. For example, it lacks modules that support high-throughput structural mutation screening, such as perturbation intensity spatial clustering, rhythm stability judgment, and reference sequence fitting difference calculation. This makes it more suitable as a preprocessing method than a terminal mutation identification module. In contrast, this invention achieves a closed-loop structure throughout the entire process, including rhythm recognition, perturbation calculation, adaptive recognition, boundary segmentation, and credibility enhancement, demonstrating technical continuity and compatibility. Summary of the Invention

[0005] The purpose of this invention is to provide a method for identifying variations in nanopore single-molecule sequencing current signals, thereby addressing some of the drawbacks and shortcomings pointed out in the background art.

[0006] The present invention addresses the aforementioned technical problems by employing the following technical solution: a method for identifying variations in nanopore single-molecule sequencing current signals, comprising: dividing the raw current signal generated by nanopore single-molecule sequencing into multiple consecutive time periods according to the order in which bases pass through the nanopore; calculating the average signal rhythm of each consecutive time period; constructing a set of standard time rhythms for standard base sequences passing through the nanopore; comparing the shift difference between the actual measured average signal rhythm and the standard time rhythm; and identifying regions with significant rhythm shift differences as potential variation regions.

[0007] Expanding outwards from the potential mutation region, extract the current signals of the neighboring regions before and after the potential mutation region, analyze the changing trends and amplitude differences of the current signals before and after the potential mutation region, quantify the disturbance difference of the current signals before and after the potential mutation region, and determine the region where the current disturbance difference has obvious directional asymmetry characteristics as the mutation region to be confirmed.

[0008] Expanding outwards from the region to be confirmed as the center, extract current signals over a wide range, calculate the disturbance intensity of the current signals at each location, analyze the variation law of the disturbance intensity within the expanded range, and determine that the region to be confirmed as the mutation has physical rationality when the disturbance intensity shows a continuous trend of gradually increasing or gradually decreasing.

[0009] For regions of variation that have physical plausibility, the theoretical current signal is determined using the corresponding unvarnished sequence. The theoretical current signal is then compared with the actual measured current signal to calculate the current difference between the two. When the current difference exceeds the normal difference range corresponding to the unvarnished sequence and cannot be reasonably explained by the unvarnished sequence, the region of variation to be confirmed is finally identified as a real variation region.

[0010] Furthermore, the method for calculating the average signal rhythm for each consecutive time period includes:

[0011] The original current signal is processed by a sliding window, and smoothed by a moving average or weighted moving average with a fixed sampling point length to obtain a smoothed current signal sequence. The smoothed current signal sequence is divided into multiple consecutive time periods. The boundary of each time period is determined according to the rate of change of the current signal, or it is divided according to the sampling interval of equal length. Each time period contains several base passage events.

[0012] For the smoothed current signal within each time period, the number of events corresponding to the current disturbance is identified. The number of events is determined by counting the number of current jump points or by counting the number of stable current intervals. Based on the total duration of each time period and the number of events within that time period, the average rhythm of each time period is calculated to characterize the average current response period per unit event or the event density per unit time.

[0013] Furthermore, the identification of the number of events includes: establishing current disturbance energy change curves for multiple consecutive time periods of the smoothed current signal, using the number of intervals between disturbance energy peaks as the basis for the number of events, and determining the number of valid events contained in each segment based on the continuity of the current disturbance.

[0014] Furthermore, the calculation process of the average rhythm includes: performing statistical analysis on the event intervals within each time period, determining whether the rhythm distribution conforms to a local normal distribution or a stable sequence pattern, and marking time periods that deviate from this pattern as abnormal segments for subsequent removal or correction.

[0015] Furthermore, the average rhythm result is used to calculate the local rhythm gradient, and based on the rate of rhythm change between adjacent time periods, the rhythm change region is marked as a candidate input for downstream variation identification analysis; the calculation of the average rhythm also includes dynamic correction of the long-term drift trend, which is achieved by establishing a global rhythm reference curve between multiple time periods and aligning the rhythm calculation results of each local segment with the reference curve.

[0016] Furthermore, the method for calculating the disturbance intensity of the current signal at each location includes:

[0017] The original current signal is preprocessed by denoising it using any one of moving average, wavelet transform, or Gaussian filtering to obtain a smooth current signal sequence. Signal data within a fixed-length window centered on each sampling point in the smooth current signal sequence is extracted. Based on the signal within the window, local mean values ​​are calculated or linear or polynomial fitting is performed to obtain local trend lines. These local mean values ​​or local trend lines are used as local references for the sampling points.

[0018] The current signal value at each sampling point is compared with the corresponding local reference to obtain the perturbation intensity of each sampling point. The perturbation intensity is the absolute difference or squared difference between the current signal value at the sampling point and its local mean or local trend line value. The perturbation intensity of all sampling points is normalized to obtain a perturbation intensity spectrum covering the entire signal sequence. Based on the normalized perturbation intensity spectrum, a perturbation intensity threshold is set, and sampling points with perturbation intensity greater than the threshold are marked as significant perturbation points.

[0019] Furthermore, for the normalized disturbance intensity spectrum, combined with the historical fluctuation patterns of the current signal, the preset threshold of the disturbance intensity is dynamically adjusted so that the threshold adapts to the noise level and fluctuation characteristics of different signal segments; the distribution of the significant disturbance points is clustered, using density clustering or continuous segmentation to merge spatially adjacent or densely distributed significant disturbance points into abnormal regions; to adapt to the noise level and background fluctuation characteristics of different signal segments, a method of dynamic threshold adjustment combined with historical fluctuation patterns is used, employing a spatial clustering mechanism to identify continuous abnormal regions; after the disturbance intensity spectrum is formed, a novel threshold adjustment formula with an exponential offset term and a variable-order logarithmic term is introduced, as follows:

[0020]

[0021] in:

[0022] For the first The adaptive disturbance intensity threshold corresponding to each sampling point; For the first The position of each sampling point in the current signal sequence; The width is half the width of the window, indicating that... The historical reference interval extends to the left and right from the center; For the perturbation intensity spectrum at position The rate of change at a point, i.e., the local first derivative; For The variance of the disturbance intensity within the centered window is used to measure noise fluctuation; The disturbance reference coefficient controls the overall level of the control threshold; This is the offset control factor, used to adjust the degree of emphasis on the near center during integration; It is a nonlinear response exponent, used to adjust the sensitivity to local noise; The basis parameter is a variable-order logarithm, used to amplify the threshold growth in high-noise regions;

[0023] By integrating the rate of change of local disturbance intensity and adjusting the distance using an exponential weighting function, the disturbance changes closer to the sampling point contribute more. Simultaneously, by introducing a variable-order logarithmic factor into the local disturbance variance, the threshold is raised in high-noise regions to avoid false alarms, while maintaining sensitivity in stable regions, thus improving recognition accuracy.

[0024] function The derivation process includes:

[0025] With sampling points Centered on the graph, we consider the changes in the perturbation intensity spectrum within its left and right windows to characterize the overall trend of the local perturbation intensity. First, we define the local variation rate of the perturbation intensity spectrum $S(x)$ as... , that is, the first derivative of the perturbation intensity at that point, represents the location and extent of the abrupt change in perturbation. To obtain... The overall trend of disturbance changes in the vicinity is analyzed, and a sliding integral mechanism is introduced, that is, in the context of... Centered on, half width is Within the window, for Integrating is performed to accumulate the total amount of disturbance change. Simultaneously, to avoid equal weighting of the contribution of all historical disturbance changes to the threshold at the current point, an exponential offset factor is introduced. It is based on Centered The type of decay function assigns higher weight to disturbance changes closer to the center, so that the threshold of the current point is mainly determined by the disturbance trend in its vicinity. This weighted integral constitutes the linear kernel of the entire formula.

[0026] However, integrating the rate of change of the perturbation alone is insufficient to address the differences in background noise levels across different signal segments; therefore, a variance factor is introduced. It means to The statistical variance of the disturbance intensity within the window, centered at [center], is used to characterize the severity of fluctuations in that segment. To avoid the threshold increasing too rapidly and nonlinearly with increasing noise, a variable-order control factor is introduced into the threshold function, i.e., [factor]. The power exponent is used as a nonlinear enhancement term for the perturbation integral result, so that the threshold is appropriately increased in the high-noise section to avoid misjudgment, while remaining sensitive in the stable section and still accurately identifying weak perturbations.

[0027] Finally, by combining the integral kernel of the disturbance change intensity and the local noise fluctuation term, the two are integrated into a power function form, while introducing an overall adjustment parameter. (Control threshold baseline level) and (Adjusting noise sensitivity) a complete adaptive disturbance intensity threshold calculation formula is constructed. The structure of this formula has a triple adjustment mechanism: ① spatial weighted integral enhances the ability to focus on continuous disturbances; ② noise variance control adapts to different fluctuation regions; ③ power exponential nonlinearity enhances the ability to respond to differences in signal characteristics, thereby achieving accurate identification and interference elimination of significant disturbance events as a whole.

[0028] Furthermore, a multi-scale analysis strategy is introduced for the calculated perturbation intensity spectrum. The perturbation intensity is calculated based on windows of different lengths, and correlation analysis is performed on the perturbation intensity spectrum at each scale to identify multi-level signal anomalies at different spatial scales.

[0029] Furthermore, based on the changing trend of the perturbation intensity spectrum over a continuous time period, change point detection or sequence segmentation methods are used to locate the start and end boundaries of the perturbation mutation, enabling the identification of the boundaries of the mutated region.

[0030] The proposed method for variant identification of nanopore single-molecule sequencing current signals forms a systematic and highly sensitive analysis process in terms of signal rhythm extraction, perturbation intensity calculation, adaptive threshold judgment, multi-scale anomaly identification, and precise boundary localization, and has the following significant advantages:

[0031] By introducing methods such as moving average and polynomial fitting to calculate local perturbation intensity, and combining this with a dynamic threshold adjustment mechanism, the system can effectively identify low signal-to-noise ratio mutation events such as single-base mutations and small insertions / deletions even in high-noise or complex backgrounds, significantly improving sequencing error control and mutation detection sensitivity. By incorporating signal history fluctuation modeling and a variable-order response function into perturbation judgment, the system effectively avoids misjudgments caused by a uniform threshold strategy, allowing the threshold to adaptively adjust according to different segments, reducing the false alarm rate in high-noise regions while maintaining high recognition accuracy in stable regions.

[0032] By constructing multi-scale perturbation spectra and conducting cross-scale consistency analysis, we can simultaneously identify local abrupt changes and domain-level slowly varying signal anomalies, providing support for the detection of structural variations (such as large segment deletions and duplications) and achieving more spatially resolving anomaly identification. Through perturbation spectrum trend modeling and employing change point detection and sequence segmentation algorithms, we can accurately define the start and end points of variation regions, solving the problems of ambiguous variation region boundaries and mis-merging of segments in traditional methods, and improving the accuracy of downstream alignment, annotation, and reliability assessment. Attached Figure Description

[0033] Figure 1 This is a simplified flowchart for identifying variations in nanopore sequencing current signals in this invention.

[0034] Figure 2 This is a diagram showing the relationship between the rhythm calculation function based on sliding window and event quantity recognition in this invention.

[0035] Figure 3 This is a flowchart of the perturbation intensity variation identification based on adaptive threshold and multi-scale analysis of the present invention.

[0036] Figure 4 This is a schematic diagram of the process for segmented rhythm recognition and variation detection of nanopore current signals in Embodiment 1 of the present invention.

[0037] Figure 5 This is a flowchart of the high-sensitivity nanopore current signal micro-perturbation detection and variation region localization in Embodiment 2 of the present invention.

[0038] Figure 6 This is a graph showing the change in energy due to current disturbance in this invention. Detailed Implementation

[0039] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0040] Combined with appendix Figure 1This invention provides a method for identifying variations in current signals during nanopore single-molecule sequencing. Based on the acquired original current signal sequence during the nanopore sequencing process, the current signal is divided into multiple consecutive time periods according to the physical order in which the base monomers in the DNA or RNA molecule pass through the nanopore. This division can be based on equally spaced sampling lengths or on the changing trend or local steady-state characteristics of the current signal to determine the boundaries of the time periods, so that each time period can effectively correspond to the passage process of a continuous set of bases. After the division is completed, the average rhythm of the current signal in each time period is calculated, where the average rhythm is used to represent the average time span corresponding to the current response triggered by a unit base in that time period. The average rhythm can be calculated by identifying the number of pore-passing events (i.e., the jump points or current patterns corresponding to the base pore-passing actions) in that time period and dividing it by the total duration of that segment. Furthermore, to obtain a comparative baseline with discriminative value, a set of standard time rhythms for standard base sequences passing through nanopore channels was constructed based on previous experimental or training data. These standard rhythms were obtained through statistical analysis of signals from multiple unvariated reference samples, representing the typical through-pore rhythm characteristics of DNA molecules from different regions passing through nanopores without sequence structural perturbation. Then, the standard time rhythms were compared with the actual measured average signal rhythms, and the time offset differences between the two were analyzed segment by segment. Specifically, the degree of local acceleration, deceleration, or rhythm imbalance of the actual rhythm relative to the standard rhythm was calculated and used as an evaluation index of rhythm consistency. When the average rhythm of the actual signal deviates significantly from the standard rhythm within a certain time period, and this deviation exceeds a preset tolerance threshold, that time period is marked as a region with significant rhythm offset and submitted as a potential candidate region for structural variation to subsequent perturbation analysis, structural verification, or sequence inference processes. This enables early identification and precise localization of regions containing sequence variations such as mutations, insertions, and deletions in nanopore current signals.

[0041] In this invention, the method for determining the tolerance threshold is as follows:

[0042] First, a standard rhythm model was constructed through statistical analysis of the unmutated reference samples to obtain the standard rhythm value corresponding to the current response induced by a unit base in each sequencing time period and its standard deviation under normal conditions. Then, based on statistical stability, a fold factor k (e.g., k=2 or k=3) is set to... The tolerance threshold for rhythm deviation is set in the form of k, where the value of k can be adjusted and optimized according to the true positive and false positive rates in the training data.

[0043] In a specific embodiment, for example, the reference rhythm of a certain segment is 4.2 ms / base, and the standard deviation is 0.3 ms / base. When the actual rhythm is detected to be 5.0 ms / base, the offset is 0.8 ms / base, which exceeds the tolerance threshold of 0.6 ms / base set with k=2. Therefore, it is marked as a region with significant rhythm offset and is further submitted to the subsequent perturbation analysis or sequence variation identification module for accurate location of potential structural variation regions. Compared with using a fixed threshold, this method is more robust and region-adaptive, and can flexibly adjust the judgment criteria according to the current response characteristics of each region, effectively improving the accuracy and recall of variation identification.

[0044] For potential variation regions identified through rhythm offset or other initial screening methods, further directional perturbation verification is performed to improve the accuracy and physical rationality of variation determination. Specifically, taking the potential variation region as the center, a neighborhood interval of a set length is extended to both sides along the front and back directions on the current signal time axis. The original or pre-processed current signal sequence within this neighborhood is extracted as the analysis object. Then, the current signal in the front neighborhood (i.e., in front of the potential variation region) and the back neighborhood (i.e., behind it) is analyzed for change trend and amplitude features. The change trend can reflect the directional perturbation characteristics of the current signal through signal slope, derivative sign change, local extreme point distribution, etc. The amplitude features can reflect the continuity and stability of signal strength through local average value, fluctuation range, standard deviation, etc. The trend direction and numerical amplitude of the signal changes in the front and back neighborhoods are compared. Further, a quantification mechanism for perturbation difference is introduced. The key current characteristic parameters of the preceding and following neighborhoods are calculated to obtain a directional perturbation difference index, which is used to measure whether the electrical signal deviation of the current central region in its preceding and following neighborhoods has obvious asymmetry. When the perturbation difference exceeds the preset asymmetry threshold, that is, when the current in the preceding and following neighborhoods shows a directional shift in trend or energy expression, it indicates that the structural perturbation of the central region is caused by actual base sequence variation (such as mutation, insertion, deletion), rather than normal signal noise or local structural oscillation. Therefore, this region is marked as a region of unconfirmed energy change variation, which serves as the key processing object for subsequent perturbation diffusion analysis, irreversible verification, or variation type identification, thereby achieving refined variation screening based on the directional difference of current perturbation.

[0045] For the variant regions to be confirmed through directional disturbance analysis, the physical rationality of the variation is further verified from the perspective of the continuity of the disturbance intensity distribution to improve the reliability and structural consistency of the variation determination. Specifically, taking the variant region to be confirmed as the center, the current signal sequence is extended to a larger range along the front and back directions to extract the current signal covering a longer time period on both sides of the region. This is used to construct the disturbance intensity sequence of the extended range. The difference between the current value of each sampling point and the mean or fitted trend value of its local window is calculated to reflect the degree of deviation of the current at that point from its background signal. Then, the disturbance intensity of all sampling points in the extended range is arranged in sequence to form a complete disturbance intensity change curve. Further continuity analysis of the curve was conducted, and the direction and gradient trend of the perturbation intensity were calculated using a sliding window to determine whether the perturbation intensity showed a gradual increasing or decreasing trend in space. When this trend persisted within a certain range without any drastic reversal or interruption, the perturbation behavior was considered to have physical continuity and reasonable diffusion. Thus, it could be determined that the region to be confirmed by the mutation was not an occasional noise or a local mutation point, but a stable perturbation center caused by base variation. The perturbation propagated forward and backward in a traceable manner, which is consistent with the physical mechanism that changes in molecular structure in nanopores have a continuous impact on current. Therefore, the region to be confirmed by the mutation was deemed to have physical reasonableness and was submitted as a high-confidence candidate region for subsequent irreversible verification or type classification steps. This further improved the identification of the mutation region from local to continuous structure confirmation, effectively enhancing the causal relationship between electrical signal anomalies and molecular structure perturbations.

[0046] For regions of variation that have been determined to be physically plausible through perturbation intensity distribution analysis, a reverse verification mechanism is further introduced. By comparing the region with the theoretical current signal corresponding to the unvariant sequence, it is verified whether the region can be reasonably interpreted by the standard sequence, so as to ultimately confirm whether it is a real variation region. Specifically, based on the reference genome or a known unvariant sequence, the base sequence information corresponding to the region of variation to be confirmed and a certain range before and after it is extracted. Based on the characteristics of the nanopore sequencing platform, the theoretical current signal sequence that the unvariant sequence should produce in this region is called or generated. This theoretical signal can be constructed by empirical current templates, statistical models or signal simulation engines, representing the standard current waveform that nanopore sequencing should output under ideal conditions where there are no sequence variations. Subsequently, the theoretical current signal is aligned point-to-point with the actual measured current signal of the region to be confirmed for mutation. Registration of the two signals on the time scale is achieved through methods such as sliding window or dynamic time comparison. Based on the registration, the current difference between the two signals is calculated. This difference can include amplitude difference, signal trend shift, cumulative residual, etc., to reflect the degree of similarity deviation between the theoretical and actual signals. Furthermore, the permissible range of current difference for unmutated signals under normal sequencing conditions is established using reference data, i.e., defining the natural error distribution interval under non-mutated conditions. When the difference between the actual current signal and the theoretical current signal significantly exceeds this normal range, and the difference lacks noise interpretability or reasonable sequencing system error, it indicates that the actual current behavior of the region cannot be reproduced or explained by the unmutated sequence. Therefore, the region to be confirmed for mutation is determined to be a true mutation region, meaning its current deviation is caused by structural perturbation due to real base variation. Finally, the region is confirmed as a mutation site or structurally abnormal region and submitted to subsequent mutation type identification, annotation analysis, or result output modules, thus achieving a closed-loop process from rhythm analysis, perturbation detection, physical verification to logical confirmation.

[0047] Combined with appendix Figure 2To address the technical aspect of calculating the average signal rhythm for each consecutive time interval, a rhythm calculation method combining sliding window smoothing and event count identification is proposed. The raw current signal sequence acquired during sequencing is processed using a sliding window. A fixed-sampling-point length moving average or weighted moving average is employed to smooth the raw current signal, reducing the impact of instantaneous noise and high-frequency fluctuations, resulting in a smoother current signal sequence that better reflects the characteristics of the base passage process. Subsequently, the smoothed current signal sequence is divided into multiple consecutive time intervals. These intervals can be divided in two ways: first, based on the rate of change of the current signal, i.e., by identifying abrupt changes in the signal slope or derivative threshold to determine the segment boundaries; second, using a fixed sampling length for equal-length division, ensuring that each time interval contains a certain number of base passage events to form an effective statistical basis. After dividing the time period, the number of events is identified for the smoothed current signal within each time period. An event refers to a significant perturbation change in the current signal representing one or more bases passing through the nanopore. The number of events can be identified in two ways: one is to count the number of jump points in the current signal, i.e., to identify the number of locations where the current value changes abruptly, reflecting the rapid transition between different states; the other is to judge based on the number of stable current intervals, i.e., to identify the number of segments in the current signal that remain stable. Each stable segment approximately corresponds to one complete base influence process. The two methods can be flexibly selected or combined according to the signal characteristics to improve the identification accuracy. After the number of events in each time period is identified, the average signal rhythm of each time period is calculated by combining the total duration information of that time period. The average rhythm can be used to characterize the average current response period of a unit event, i.e., to represent the average passage time corresponding to each base or structural event in that time period. It can also be used to inversely represent the event density by calculating the number of events contained in a unit of time. This rhythm parameter serves as an important basic quantitative indicator for judging rhythm shifts, rhythm inconsistencies, and regional anomalies in subsequent variation identification, which helps to achieve preliminary identification of potential variation regions and early warning of structural perturbations.

[0048] Combined with appendix Figure 6As shown, to improve the stability and accuracy of event recognition, an event quantity recognition method based on disturbance energy change curves is proposed. This method processes a pre-processed smooth current signal sequence, dividing the entire signal sequence into multiple continuous time periods. For each time period, a disturbance energy change curve is constructed for that current signal. The disturbance energy change curve is used to characterize the fluctuation intensity distribution of the current signal on the time axis. A fixed-length symmetrical window is extracted centered on each sampling point, and the local mean is calculated from the current values ​​of all sampling points within the window as a reference baseline. Subsequently, the deviation of the current at each point within the window from the reference value is calculated, and the deviation values ​​are squared. Finally, all squared deviation values ​​are summed to form the disturbance energy value corresponding to the current sampling point. This disturbance energy value is arranged along the time axis to form a complete disturbance energy change curve, which is used to identify current disturbance events.

[0049] This method reflects the spatial evolution of the current perturbation intensity caused by the passage of a local structure through a nanopore. Further analysis of the energy peak distribution characteristics based on this perturbation energy curve is conducted. By setting reasonable minimum peak spacing, minimum peak amplitude difference thresholds, and response widths, effective local peak points on the energy curve are identified. These peak points represent strong response events during the current perturbation process. The interval between two adjacent peaks can be considered as the time interval for an independent base or structural unit to pass through the nanopore. The number of intervals between peaks is then counted, providing an estimate of the number of events within the current time interval. Unlike traditional methods that rely solely on current jump points, this approach analyzes the continuity and morphological characteristics of the perturbation energy, maintaining event identification stability while exhibiting stronger noise resistance and adaptability to complex signal structures. Finally, the average rhythm or rhythm density is calculated based on the number of events identified in each segment and the segment length, providing reliable basic data for subsequent rhythm shift analysis of variation regions, structural consistency assessment, and physical rationality judgment.

[0050] To further improve the accuracy and robustness of rhythm analysis, a rhythm distribution regularity judgment mechanism is introduced into the calculation of average rhythm. Specifically, after identifying the number of events within each time period, the interval sequence between adjacent events within that time period is obtained. The event interval reflects the temporal distribution characteristics of current signal perturbations during the passage of continuous bases or structural units through nanopores. Subsequently, statistical analysis is performed on this event interval sequence, calculating basic statistical characteristic parameters such as mean, variance, skewness, and kurtosis to measure its central tendency and fluctuation characteristics. Then, the interval distribution is further fitted or tested to determine whether it conforms to the characteristics of a local normal distribution, i.e., whether the event intervals exhibit a stable symmetrical distribution structure around a certain average period, or whether it satisfies the stable sequence regularity, i.e., the variation of event intervals on the time axis. The method assesses whether the sequence exhibits low fluctuations, no abrupt changes, and near-stationary characteristics. When the event intervals within a certain time period show abnormal fluctuations, such as a significant increase in variance, severe deviation of skewness from zero, frequent occurrences of excessively short or long intervals, or significant discrepancies after a normality test, the rhythm distribution within that time period is determined to be abnormal. This indicates that the segment is affected by factors such as local strong noise, unstructured disturbances, or signal drift. Therefore, it is marked as an abnormal segment and enters the subsequent elimination or correction process. Specific processing methods may include excluding the time period from subsequent variation analysis, interpolating and smoothing or replacing its rhythm values, or re-dividing the time period for corrective calculation. This method improves the reliability of average rhythm estimation by analyzing the regularity of the event interval distribution structure and effectively avoids abnormal segments interfering with the accuracy of downstream variation location.

[0051] To enhance the sensitivity of average rhythm to local anomalies and its robustness to systematic drift in identifying real variations, a dual mechanism based on rhythm gradient analysis and global reference alignment is proposed. After obtaining the average rhythm results for multiple consecutive time periods, the rhythm difference between adjacent time periods is calculated, and this difference is divided by the time or position interval between the two periods to obtain the local rhythm gradient, which is used to quantify the rate and amplitude of rhythm changes. When the rhythm gradient of a certain region is significantly higher than that of adjacent regions, it indicates a sudden change in rhythm at that point, corresponding to a change in local sequence structure or the influence of abnormal bases. Therefore, this rhythm-changed region is marked as a candidate variant region and submitted to the downstream perturbation intensity analysis or structural rationality judgment module as input, thereby realizing early variant clue identification based on rhythm change trends. Secondly, to overcome the problem of overall rhythm baseline shift caused by factors such as changes in channel electrochemical environment, sequencing speed fluctuations, or equipment drift during signal sequencing, a global rhythm reference alignment mechanism is further introduced. That is, based on the average rhythm results of all time periods, a rhythm baseline curve is constructed. The baseline curve can be obtained by multi-segment moving average, curve fitting, or spline smoothing, and is used to represent the overall trend of rhythm changes in the entire signal sequence. Then, the average rhythm results of each local time period are compared with the baseline curve, and time periods that deviate significantly from the global trend are corrected or their weights adjusted to ensure that the rhythm calculation results have relatively stable global consistency. By combining this method of local sensitive mutation detection with overall baseline correction, the sensitivity of variant region identification is improved while effectively controlling rhythm misjudgment caused by systematic drift.

[0052] Combined with appendix Figure 3To accurately identify abnormal response characteristics in current signals caused by sequence variations or structural disturbances, a disturbance intensity extraction method based on local reference difference calculation is proposed to achieve quantitative characterization of the disturbance degree at each location in the signal. The original current signal is preprocessed using moving average, wavelet transform, or Gaussian filtering to eliminate interference from transient changes, instrument background noise, or short-term high-frequency disturbances, resulting in a smoothed current signal sequence. Then, a window of data of a certain length is extracted centered on each sampling point in the smoothed current signal sequence, serving as the local analysis range for that point. Statistical modeling of the signal data within this window is performed, calculating the local mean or constructing a local trend line through linear regression, polynomial fitting, etc. The local mean and trend line reflect the expected current level of that point within its neighborhood, serving as local reference signals. Finally, the current value at the current sampling point is compared with the corresponding local reference value, and the difference is calculated. This difference represents the disturbance intensity at that point, which can be determined by the difference between the current value and the local mean or fitted value. The absolute difference or squared difference is calculated to reflect the degree of deviation of the point relative to the local background. After obtaining the perturbation intensity of all sampling points, the perturbation intensity values ​​of the entire sequence are further normalized to make the perturbation intensity of different regions have uniform dimensions and comparability. The normalization method can be mean-variance standardization or maximum-minimum linear scaling, etc., to finally form a perturbation intensity spectrum covering the entire signal sequence, which serves as the basic indicator for anomaly point screening. Further, a perturbation intensity threshold is set based on the normalized perturbation intensity spectrum. This threshold can be a globally fixed value or dynamically adjusted according to the local noise level of the signal to determine whether the current deviation is significant. When the perturbation intensity of a certain sampling point exceeds the threshold, the point is considered a significant perturbation point, indicating that it is affected by structural variation, sequence change or other abnormal interference, and can be used as a key signal feature input for subsequent candidate region detection.

[0053] To further improve the accuracy and adaptability of current signal disturbance intensity identification, after constructing and normalizing the disturbance intensity spectrum of the original current signal, a technical solution combining historical fluctuation patterns for dynamic threshold adjustment is proposed. This solution analyzes the local fluctuation characteristics of the disturbance intensity spectrum in different regions and adaptively adjusts the intensity threshold used to determine significant disturbance points. This ensures that the threshold can respond sensitively to noise backgrounds and disturbance behaviors in different signal segments, avoiding false alarms in high-noise areas or missing true variation signals in low-disturbance areas. Specifically, after the disturbance intensity spectrum is formed, the threshold is used as an arbitrary value in the signal sequence... Centered on a sampling point, a window width of [value] is set with that point as the center. Within a reference interval, the rate of change of the disturbance intensity spectrum (i.e., the local first derivative) is calculated. An exponential decay function is introduced to weight the spatial contribution of the rate of change integral, making the disturbance change near the center point contribute more to the threshold, thus forming the weighted disturbance integral expression for that point. Simultaneously, the variance of the disturbance intensity value is calculated within this local window as an indicator of the fluctuation level of the signal segment at that point. The result of the rate of change integral and the variance result are then combined into a dynamic threshold formula with a nonlinear control term, as follows:

[0054]

[0055] in, For the first Dynamic threshold for each sampling point This represents the position of the sampling point in the signal sequence. The width of the window represents the range of the historical reference interval, which is expanded horizontally. Indicates the perturbation intensity spectrum at position The rate of change at a given point, i.e., the local first derivative at that location, is used to reflect the degree of abrupt change in the perturbation trend. The variance of the disturbance intensity within the central window of the sampling points is used to characterize the level of local noise fluctuation. The disturbance baseline control coefficient determines the overall threshold level. This is the exponential offset factor, which controls the weighting priority of perturbations near the center point in the integral. To respond to the exponent, the nonlinear relationship between perturbation changes and threshold increases is adjusted. The variable-order logarithmic basis parameter is used to amplify the threshold rise in high-noise regions, thus forming an adaptive threshold that incorporates both the influence of local rate of change and sensitive local noise control. Through this dynamic threshold mechanism, a precise judgment baseline is set for each sampling point in the normalized perturbation intensity spectrum. Based on this, sampling points with perturbation intensities greater than the corresponding threshold are marked as significant perturbation points. Density clustering or continuous segmentation methods are then used to analyze the spatial distribution of significant perturbation points. When multiple significant perturbation points exhibit spatial proximity or high concentration on the time axis, they are merged into a single anomalous region. This achieves adaptive discrimination of significant perturbation points while further constructing candidate regions of variation with structural coherence and physical interpretability, improving the accuracy, stability, and robustness of the overall detection system in complex sequencing signal backgrounds.

[0056] To more comprehensively capture current perturbation behavior at different scales and thus improve the ability to identify complex variable signals, a multi-scale analysis strategy based on perturbation intensity spectrum is proposed. Specifically, after denoising and smoothing the original current signal and initially constructing the perturbation intensity spectrum, a perturbation intensity calculation process with multiple window lengths is introduced. Local perturbation analysis is performed on each sampling point in the signal sequence using short, medium, and long window lengths, respectively. That is, at each scale, neighborhood current data is extracted centered on its respective window, and the perturbation intensity value at the corresponding scale is calculated based on local mean difference or trend fitting residuals. This yields multiple sets of perturbation intensity spectra with different scales and resolutions on the same current signal sequence. Shorter scales are more sensitive to small, rapid signal mutations and are suitable for capturing single-base mutations or short insertion / deletion variants, while longer scales are suitable for identifying slowly changing signal trends or large-segment structural anomalies, exhibiting stronger background denoising capabilities and variant morphology integration characteristics. Further, multi-scale perturbation... After the spectrum is generated, correlation analysis is performed on the correspondence between the spectra at different scales. By sliding comparison or position alignment, it is evaluated whether the changes in the perturbation intensity of the same position or adjacent regions at different scales have a consistent trend or consistency in strength. When a position shows a synchronous increase or synchronous abrupt change in perturbation intensity at multiple scales, and its local perturbation behavior shows positive correlation or cooperative change characteristics at each scale, it can be considered that the position has consistent perturbation characteristics at multiple scales, which further enhances its credibility as a candidate region for mutation. Conversely, if a position shows anomalies only at a single scale and remains stable at other scales, it is a false positive caused by random noise or algorithm edge effects, which can be effectively eliminated through a multi-scale information mutual verification mechanism. This method introduces multi-level feature integration at spatial scales, so that mutation identification not only relies on the judgment of a single resolution, but also on cross-scale interference verification and signal consistency assessment, which improves the detection capability of base variation, structural mutation and signal distortion in complex and physiological backgrounds.

[0057] To achieve precise localization of the boundary of the variation region, a boundary identification mechanism based on the changing trend of the perturbation intensity spectrum is proposed. Specifically, after constructing and normalizing the perturbation intensity spectrum, trend analysis is performed on the perturbation intensity spectrum corresponding to continuous time periods in the current signal using change point detection or sequence segmentation methods. By sliding along the time axis to analyze the local mean, variance, slope, and other statistical characteristics of the perturbation intensity, the time points where significant abrupt changes in the spectrum value occur within a local range are identified. These abrupt changes often represent the beginning or end of the perturbation, i.e., the physical boundary of the potential variation region. To improve the sensitivity and robustness of change point identification, a double-sided sliding window comparison strategy is introduced, that is, for each candidate position on the perturbation spectrum sequence, its left window is compared with its right window. The statistical differences in the window indicators are analyzed. If the difference exceeds a set threshold, the location is determined to be a change point. Simultaneously, a dynamic segmentation method can be used to optimally segment the entire perturbation intensity spectrum, ensuring that the signal statistical characteristics within each sub-interval remain relatively stable, while the characteristic differences between adjacent sub-intervals are maximized. This achieves the segmentation and boundary manifestation of the perturbation structure. Combined with the above change point detection results, the start and end points of the perturbation mutation can be effectively marked, thereby clarifying the boundary range of the mutation region. This provides an accurate signal basis for subsequent structure type determination, mutation reliability assessment, and base-level localization, realizing a closed-loop transition from coarse-grained mutation identification to fine-grained boundary labeling, and improving the accuracy, boundary integrity, and physical consistency of the mutation identification process in nanopore sequencing current signals.

[0058] Example 1:

[0059] Combined with appendix Figure 4In this embodiment, the raw data of nanopore sequencing of a gene fragment contains a continuous sequence of sampling points in a sequencing channel. The total length is 30,000 sampling points, the sampling frequency is 1,000 times per second, and the total sequencing time is 30 seconds. The raw current signal contains typical base perturbation behavior. The average duration of the current fluctuation caused by each base passing through the nanopore is about 4 milliseconds. Theoretically, one base event spans about 4 sampling points. First, the original current signal was processed using a sliding window method. To improve signal quality and reduce sudden noise interference, a sliding average method with a length of 5 sampling points was used for smoothing. That is, the value of each point is the average of itself and the two points before and after it. After smoothing, the signal still retains obvious perturbation patterns in the base through-hole segment, while the noise spikes are significantly reduced. Then, the smoothed signal was divided into multiple consecutive time periods using a fixed equal-length strategy, with each time period consisting of 1,000 sampling points (corresponding to 1 second), resulting in 30 segments. Each segment theoretically contains approximately 250 base events. Next, the number of events in each smoothed signal segment was identified. In the first time period, the current jump point identification method was used. After the signal was first-order differencing, the number of points with a difference value greater than a certain jump threshold (e.g., 0.3 nA) was counted as the perturbation boundary. If there is a certain interval between consecutive jump points (e.g., more than 4 points), they can be identified as two independent events. In the first segment, 220 jump points met the condition, indicating that there were 220 valid bases in this segment. The total duration of this segment is 1 second, so the average rhythm can be calculated as 1 second / 220 ≈ 4.55 milliseconds / event, representing the current response period per unit event. In the second segment, if 250 events are identified, the average rhythm is 4.00 milliseconds / event; in the third segment, 270 events are identified, so the average rhythm is 3.70 milliseconds / event, forming the following rhythm sequence: 4.55, 4.00, 3.70... This continuous change in rhythm can be used to assess whether there is a rhythm shift signal. For example, if the rhythm suddenly drops to 2.80 milliseconds / event in the fifth segment, significantly deviating from the rhythm trend of the previous segments, it suggests that there is a structural variation or base insertion in this segment. In addition, the event density per unit time can also be calculated. For example, the event density is 220 times / second in the first segment and 270 times / second in the third segment. The increase in rhythm density means that the local sequence structure causes the through-hole rhythm to be abnormally accelerated. This information can be compared with the reference rhythm of the standard unvariant sample to further assist in identifying variation clues.

[0060] For the current signal that has been smoothed and divided into 30 time periods, segments 6 to 10 were selected for further analysis. Segment 6 corresponds to a sampling point range of 6000 to 7000. The number of events was identified using a perturbation energy curve method. Specifically, a perturbation energy curve was constructed for each time period. This curve was obtained by calculating the squared difference of each sampling point relative to its local mean, and integrating the perturbation intensity within each time period using a sliding window to obtain a continuous energy intensity trajectory. In the perturbation energy curve of segment 6, 11 obvious local peaks were observed, with the interval between adjacent peaks approximately 40 sampling points, corresponding to about 4 base events. This indicates that segment 6 contains approximately 44 base events. Statistical peak... The number of intervals, i.e., the number of events, was 44. Entering segment 7, the peak density of disturbance energy increased significantly, identifying 62 valid peaks, and the number of events increased significantly. Subsequently, the intervals between events in each segment were statistically analyzed. In segment 6, the mean interval was 22.7 sampling points, the standard deviation was 2.4, the skewness was close to 0.1, the distribution curve was approximately bell-shaped, and the goodness of fit was better than 0.93, indicating that it conformed to a local normal distribution and belonged to a stable rhythm segment. However, in segment 8, the standard deviation of the intervals widened to 6.7, the skewness reached 1.2, and extreme values ​​of less than 10 and greater than 35 frequently appeared. Furthermore, it failed the normality test through the Shapiro-Wilk test, and the goodness of fit was less than 0.65. Therefore, segment 8 was marked as rhythmic anomaly. The first segment provides a basis for subsequent elimination or rhythm correction. Next, the local rhythm gradient between segments 6 and 10 is calculated. The average rhythm of segment 6 is 22.7 sampling points / event, segment 7 is 17.5, segment 8 is 15.2, segment 9 rises to 20.8, and segment 10 stabilizes at 21.5. First-order differencing yields rhythm gradients of -5.2, -2.3, +5.6, and +0.7, respectively. A significant negative gradient is observed between segments 6 and 8, indicating the highest rhythm mutation rate and suggesting a high-risk region for local mutations. Therefore, segments 7 and 8 are uniformly labeled as rhythm mutation regions and used as candidate mutation inputs. Furthermore, to avoid overall rhythm shifts due to system drift or changes in the micro-conditions of the sequencing channels... To correct misjudgments, the average rhythm value of each of the 30 segments was extracted to construct a rhythm reference curve. A cubic spline function was used for smooth fitting to form a global baseline curve. The difference between the average rhythm of each segment and the reference curve was used as an offset for correction. In this example, the original average rhythm of segment 8 was 15.2, which was corrected to 17.1 after global baseline adjustment, making it closer to the rhythm level of the preceding and following segments and avoiding misjudgment as a system mutation. Finally, by shifting the event identification method from jump points to disturbance energy peaks, rhythm stability and normality analysis were used to identify rhythm abnormal segments. Combined with the dual correction mechanism of local gradient mutation and global trend regression, segments 7-8 were accurately marked as significantly offset rhythm abnormal regions and identified as candidate mutation regions.

[0061] Example 2:

[0062] Combined with appendix Figure 5 The key to highly sensitive identification of weak but stable structural perturbations in current signals lies in the calculation of perturbation intensity at the sampling point level and the labeling of significant perturbation points. A raw nanopore sequencing current signal sequence with 20,000 sampling points was obtained, corresponding to a sampling frequency of 4,000 times per second and a total duration of 5 seconds. The signal amplitude ranged from 60 pA to 100 pA. In the raw signal, due to interference from factors such as molecular pore processes, liquid oscillations, and instrument background noise, the signal exhibited relatively violent short-term fluctuations and irregular spikes, making it difficult to directly identify weak structural perturbations. Therefore, the raw signal was first preprocessed. Gaussian filtering was selected to preserve the low-frequency structural features of the signal, and convolution was performed using a 5-point wide Gaussian kernel function with a standard deviation of 2. The smoothed current signal sequence after filtering effectively suppressed high-frequency perturbations while retaining the mid-to-low-frequency variation characteristics representing the influence of structural pores.

[0063] Then, taking each sampling point in the smoothed signal sequence as the center, a local window of 21 points is extracted, consisting of 10 points to the left and 10 points to the right, serving as the local analysis interval. Within this local interval, two methods are used to calculate the local reference: one is the local mean, which directly calculates the average of the 21 points as the local expected current value for the current sampling point; the other is polynomial fitting, where a second-order polynomial is used to fit the signal within the window to obtain the trend prediction value for the current point. Taking sampling point number 7,000 as an example, its local window current data is as follows: The mean current is 74.7 pA, and the fitted value is 74.5 pA. If the current at this point is 77.6 pA, then the disturbance intensity (in terms of the difference of squares) is (77.6-74.7). 2 =8.41pA 2 Alternatively, the absolute difference method yields 2.9 pA.

[0064] After performing the same processing on all sampling points in the entire sequence, a complete perturbation intensity sequence is formed. Since the background fluctuations differ across segments, the perturbation intensity spectrum is normalized to eliminate the influence of absolute dimensions. Here, the z-score normalization method is used, which standardizes each value using the mean and standard deviation of the entire perturbation intensity. Statistical analysis yields the mean perturbation intensity of the entire sequence as follows: The standard deviation is Then, the normalized perturbation intensity at point 7,000 is (8.41-1.45) / 0.88≈7.92, indicating that it is much higher than the average perturbation background.

[0065] Then, a disturbance intensity threshold is set, using an empirical value. In principle, a threshold is set as follows: All normalized perturbation intensity spectra were compared point by point, and all sampling points with values ​​greater than the threshold were marked as significant perturbation points. In this example, approximately 420 significant perturbation points were marked, accounting for about 2.1%. These points were mainly concentrated in the structural switching region where current fluctuations were abnormally violent.

[0066] Furthermore, by analyzing the spatial distribution of these significant disturbance points, it can be found that they often do not appear in isolation, but rather in clusters. For example, there are more than 20 high disturbance points in the range of sampling points 6,900 to 7,150. Combined with the upstream and downstream rhythm shift characteristics, these points can be further identified as potential structural variation areas.

[0067] Furthermore, a dynamic threshold adjustment mechanism for the perturbation intensity spectrum and a clustering reduction method for significant perturbation points are introduced to achieve highly robust identification of variant regions. Previously, a normalized perturbation intensity spectrum was established using a sliding window and polynomial fitting, and a large number of preliminary significant perturbation points were labeled.

[0068] In this example, the region from sampling points 6,500 to 7,500 in the preceding signal sequence is selected, with a length of 1,001 points. For each sampling point within this region... Constructing an adaptive threshold Use the following adjustment formula:

[0069]

[0070] The variables are defined as follows, and the following parameters are used as example calculation inputs:

[0071] Setting the window half-width to 20 indicates a total window length of 41 points;

[0072] : is the disturbance reference control coefficient, which is taken as in this example. ;

[0073] : Offset control factor, controls the weighting degree of the center point, set ;

[0074] : Nonlinear response exponent, taken in this example ;

[0075] : Variable-order logarithmic basis, let ;

[0076] The first-order difference of the perturbation intensity spectrum can be approximated by the five-point central difference method.

[0077] :by The variance of the local disturbance intensity centered on the sphere reflects the level of local fluctuations.

[0078] Taking the 7,000th sampling point as an example, 20 points before and after it are extracted to form a window. The variance of the disturbance intensity is calculated from this window. The integral part of the rate of change of the disturbance intensity spectrum is obtained by numerical calculation as follows:

[0079]

[0080] Then its variable-order logarithmic factor is The final exponent term is Therefore, the adaptive threshold at this point is:

[0081]

[0082] If the normalized value of the disturbance intensity at a point is 4.1, significantly higher than the threshold of 3.48, it is marked as a significant disturbance point; conversely, if the integral of the rate of change in a neighboring region, such as point 7,005, is lower (4.8) and the variance is 0.18, the threshold will become:

[0083]

[0084] If the actual perturbation is 2.95, it is below the threshold and is not marked as abnormal.

[0085] After calculating the dynamic threshold for all sampling points, density clustering (such as DBSCAN) is further used to aggregate significant disturbance points. With a minimum cluster size of 8 points and a maximum distance of 15 sampling points, the clustering algorithm aggregates a continuous and dense cluster of disturbance points between 6,920 and 7,040, identifying it as anomalous region 1. Another anomalous cluster is identified between 7,180 and 7,210. Regions with scattered or widely spaced edge points (such as single-point or double-point sudden disturbances) are not included in the cluster centers by the clustering algorithm and are considered background fluctuation pseudo-signals and excluded.

[0086] In the sampling interval from 6,500 to 7,500, after completing the clustering of significant perturbation points and abnormal regions based on dynamic threshold identification, a multi-scale perturbation intensity spectrum construction strategy was introduced to improve the bidirectional identification capability of local small-amplitude abrupt changes and cross-domain slowly varying signals. At the same time, a change point detection mechanism was adopted to perform high-precision segmentation of the physical boundaries of the clustered regions, so that the entire identification process has good structural coherence and spatial resolution.

[0087] First, based on the 6,500 to 7,500 segment, three sets of perturbation intensity spectra at different scales were constructed, corresponding to the window lengths respectively. , , Each sampling point corresponds to an analysis granularity of small-scale (base level), medium-scale (local base sequence), and large-scale (domain level). The specific method is as follows: Using each sampling point as the center, extract... , , A window of length is used, and within each window, the perturbation intensity at the current sampling point is calculated based on the local mean or second-order polynomial trend, resulting in a complete perturbation spectrum sequence at each scale. Taking sampling point 7000 as an example, its original current is 77.6 pA, and the local mean is 75.1 at the small scale, 74.6 at the medium scale, and 74.3 at the large scale. The corresponding perturbation intensities are 2.5, 3.0, and 3.3 pA, respectively. After normalization, these values ​​are 2.9, 3.2, and 3.6 pA, indicating that as the scale increases, the assessment of the local deviation of the perturbation at this point gradually strengthens.

[0088] Subsequently, a correlation analysis was performed on the three-scale spectra to align their positions. Specifically, based on the position of each sampling point, a vector was formed by taking the perturbation intensity values ​​at each of the three scales. The trends of these vectors on the signal sequence were calculated by sliding. Within the window of 7000 to 7020, the multi-scale spectrum showed a consistent upward trend at points 7006 and 7014, with a Pearson correlation coefficient as high as 0.94 in this small window, suggesting a consistent perturbation enhancement behavior at these positions. However, in the interval of 7035 to 7050, the large-scale spectrum increased while the small-scale spectrum was relatively flat, with a correlation coefficient of only 0.42, indicating a slow background drift and atypical variation behavior.

[0089] Based on the above results, we further focused on the abnormally high segment between 7000 and 7030, and used a change point detection method to locate its boundary. In this example, we used Cumulative Sum Change Detection (CUSUM) combined with the sliding variance mutation method, with a sliding window of 5 points, to monitor the rate of change of the mean and variance in the perturbation intensity spectrum. Two abrupt jumps in the statistical characteristics of the perturbation were detected at 7003 and 7026. The first was a jump in the perturbation mean from 1.9 to 3.4, and the second was a drop in the mean from 3.6 back to 2.1, with the variance decreasing significantly at the same time. Therefore, this segment was marked as the perturbation mutation region, with its physical boundary clearly defined as the sampling points from 7003 to 7026, a length of 24 sampling points, corresponding to a time of 6 milliseconds, theoretically covering approximately 6 base events.

[0090] Thus far, in the 6,500 to 7,500 segment, not only were composite perturbation features with both local mutations and macroscopic background shifts successfully extracted using multi-scale windows, but the start and end positions of typical variation regions (e.g., 7003–7026) were also accurately located using change point detection methods, laying a spatial foundation for subsequent candidate sequence comparison, theoretical current fitting, and structural confirmation.

[0091] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A method for identifying variations in nanopore single-molecule sequencing current signals, characterized in that... include: The raw current signal generated by nanopore single-molecule sequencing is divided into multiple consecutive time periods according to the order in which the bases pass through the nanopore. The average signal rhythm of each consecutive time period is calculated, and a set of standard time rhythms for standard base sequences passing through the nanopore are constructed. The shift difference between the actual measured average signal rhythm and the standard time rhythm is compared, and the rhythm shift difference region is identified as a potential variant region. Expanding outwards from the potential variation region, extract the current signals of the neighboring regions before and after the potential variation region, analyze the changing trends and amplitude differences of the current signals before and after the region, quantify the disturbance difference of the current signals before and after the region, and determine the region where the current disturbance difference has directional asymmetry characteristics as the variation to be confirmed region. Expanding outwards from the region to be confirmed as the center, extract current signals over a wide range, calculate the disturbance intensity of the current signals at each location, analyze the variation law of the disturbance intensity within the expanded range, and determine that the region to be confirmed as the mutation has physical rationality when the disturbance intensity shows a continuous trend of gradually increasing or gradually decreasing. For regions of variation that have physical plausibility, the theoretical current signal is determined using the corresponding unvarnished sequence. The theoretical current signal is then compared with the actual measured current signal to calculate the current difference between the two. When the current difference exceeds the normal difference range corresponding to the unvarnished sequence and cannot be reasonably explained by the unvarnished sequence, the region of variation to be confirmed is finally identified as a real variation region.

2. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 1, characterized in that... The method for calculating the average signal rhythm for each consecutive time period includes: The original current signal is processed by a sliding window, and smoothed by a moving average or weighted moving average with a fixed sampling point length to obtain a smoothed current signal sequence. The smoothed current signal sequence is divided into multiple consecutive time periods. The boundary of each time period is determined according to the rate of change of the current signal, or it is divided according to the sampling interval of equal length. Each time period contains several base passage events. For the smoothed current signal within each time period, the number of events corresponding to the current disturbance is identified. The number of events is determined by counting the number of current jump points or by counting the number of stable current intervals. Based on the total duration of each time period and the number of events within that time period, the average rhythm of each time period is calculated to characterize the average current response period per unit event or the event density per unit time.

3. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 2, characterized in that... The identification of the number of events includes: establishing current disturbance energy change curves for multiple consecutive time periods of the smoothed current signal, using the number of intervals between disturbance energy peaks as the basis for the number of events, and determining the number of valid events contained in each segment based on the continuity of the current disturbance.

4. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 3, characterized in that... The calculation process of the average rhythm includes: statistically analyzing the event intervals within each time period, determining whether the rhythm distribution conforms to a local normal distribution or a stable sequence pattern, and marking time periods that deviate from the pattern as abnormal segments for subsequent removal or correction.

5. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 4, characterized in that... The average rhythm result is used to calculate the local rhythm gradient, and the rhythm change rate between adjacent time periods is used as the basis to mark the rhythm change region as the candidate input for downstream variation identification analysis. The calculation of the average rhythm also includes dynamic correction of the long-term drift trend. The correction is achieved by establishing a global rhythm reference curve between multiple time periods and aligning the rhythm calculation results of each local segment with the reference curve.

6. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 1, characterized in that... The method for calculating the disturbance intensity of the current signal at each location includes: The original current signal is preprocessed by denoising it using any one of moving average, wavelet transform, or Gaussian filtering to obtain a smooth current signal sequence. Signal data within a fixed-length window centered on each sampling point in the smooth current signal sequence is extracted. Based on the signal within the window, local mean values ​​are calculated or linear or polynomial fitting is performed to obtain local trend lines. These local mean values ​​or local trend lines are used as local references for the sampling points. The current signal value at each sampling point is compared with the corresponding local reference to obtain the perturbation intensity of each sampling point. The perturbation intensity is the absolute difference or squared difference between the current signal value at the sampling point and its local mean or local trend line value. The perturbation intensity of all sampling points is normalized to obtain a perturbation intensity spectrum covering the entire signal sequence. Based on the normalized perturbation intensity spectrum, a perturbation intensity threshold is set, and sampling points with perturbation intensity greater than the threshold are marked as significant perturbation points.

7. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 6, characterized in that... The normalized disturbance intensity spectrum is combined with the historical fluctuation pattern of the current signal to dynamically adjust the preset threshold of the disturbance intensity, so that the threshold adapts to the noise level and fluctuation characteristics of different signal segments. The distribution of the significant disturbance points is clustered using density clustering or continuous segmentation methods to merge spatially adjacent or densely distributed significant disturbance points into anomaly regions.

8. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 7, characterized in that... For the calculated perturbation intensity spectrum, a multi-scale analysis strategy is introduced to calculate the perturbation intensity based on windows of different lengths, and correlation analysis is performed on the perturbation intensity spectrum at each scale to identify multi-level signal anomalies at different spatial scales.

9. The method for identifying variations in nanopore single-molecule sequencing current signals according to claim 8, characterized in that... To identify the changing trend of the perturbation intensity spectrum over a continuous time period, change point detection or sequence segmentation methods are used to locate the start and end boundaries of the perturbation abrupt change, thereby enabling the identification of the boundaries of the variability region.

Images