A method and medium for analyzing sparse microbiome time series data

Abstract

This invention discloses a sparsity-based microbiome time-series data analysis method and medium, relating to the field of microbial analysis technology. The method includes acquiring a microbiome time-series data matrix and initializing a centroid set and a non-negative time dimension weight vector. Subsequently, iterative optimization is performed with a minimum objective function as the constraint: samples are assigned to the nearest centroid to generate cluster partitions; local and global medians are extracted to generate inter-cluster sums of squares for each time dimension; a soft threshold operator is used to apply a regularization penalty to the inter-cluster sums of squares, redundant time dimension weights are truncated to zero and normalized updates are performed; finally, centroids are redefined until convergence. This invention overcomes the sensitivity of traditional mean-based evaluation systems to outliers, constructing a highly robust sparsity mechanism in a sparse, high-noise biological space, and accurately locking the true evolutionary time window.

Description

Technical Field

[0001] This application belongs to the field of microbial analysis, specifically relating to a method and medium for analyzing sparse microbiome time-series data. Background Technology With the rapid popularization of high-throughput sequencing technology, the mining and analysis of longitudinal microbiome data has become a core technical means to reveal the evolutionary patterns of host health and the pathogenesis of diseases. During long-term continuous sampling, researchers typically need to perform cluster analysis on microbial time-series data to accurately identify microbial groups with co-evolutionary characteristics. However, in the lengthy time-series acquisition period, not all physical sampling points make a substantial contribution to the evolutionary classification of microbial communities; a large number of time dimensions are actually filled with statistically insignificant background fluctuations. To suppress the interference of these redundant time dimensions and reduce data dimensionality, the industry commonly introduces sparse clustering algorithms based on regularization penalty mechanisms into traditional clustering frameworks, attempting to achieve automatic selection of key time nodes through adaptive scaling and truncation of feature weights.

[0002] In existing sparse clustering algorithms, a feature importance evaluation metric, the Between-Cluster Sum of Squares (BCSS), is typically constructed. A soft threshold penalty is then applied based on this metric to update the weights of each dimension. This traditional architecture fundamentally relies heavily on calculating the arithmetic mean of the dataset to establish the absolute center of the clusters and the global reference point. However, in real-world sequencing environments, microbiome data inherently exhibits extreme sparsity, non-normal distribution, and extreme heavy-tailed effects caused by high-frequency sequencing artifacts. When sporadic, extremely high-abundance outliers exist in the data space, the arithmetic mean in the underlying calculation is severely shifted due to strong gravitational pull, causing the BCSS-based calculation to fail to objectively reflect the true cluster dispersion of the data over a specific time dimension.

[0003] The distortions in existing technologies when dealing with noisy biological data directly lead to the failure of the sparsity penalty mechanism. Specifically, the inter-cluster sum of squares distorted by outliers, after being penalized by the regularization operator, will output completely erroneous feature weight gradients to the algorithm. This not only makes the model highly susceptible to overfitting and retains invalid time points containing a large amount of noise, but also introduces false negative errors, treating key time windows that truly signify gradual biological differentiation as redundant features and forcibly eliminating them. Summary of the Invention

[0004] Firstly, in view of the shortcomings of the existing technology, the purpose of this application is to provide a sparse microbiome time series data analysis method, which solves the technical problems raised in the background technology.

[0005] The objective of this application can be achieved through the following technical solutions: A method for analyzing sparse microbiome time series data, the method comprising: Obtaining the time-series data matrix of the microbiome The time series data matrix Include Samples of each operational classification unit in The values ​​at each aligned time point Characterizing the time series data matrix The first in Samples of each operational classification unit; Get the set number of clusters and regularization parameters ; In response to the acquired time series data matrix and the number of clusters From the time series data matrix Extract initialization from Each center point set of center points And based on the time series data matrix The time dimension attribute initialization includes Each weight value Non-negative time dimension weight vector ; In response to the acquired time series data matrix The initial set of center points and the non-negative time dimension weight vector The algorithm iteratively solves the problem by minimizing the following objective function, determining whether the parameter changes satisfy a preset convergence condition. If the result indicates that the preset convergence condition is met, the current clustering partition is output. The set of center points and the non-negative time dimension weight vector As the optimal combination of solutions to the objective function; The iterative solution includes: Obtain the samples of the operation classification unit To each of the aforementioned center points The weighted distance is used to classify the samples of each of the operation classification units. Uniquely assigned to the nearest center point The corresponding cluster In the process, the clustering partition is generated. ; In response to the generated cluster partition From the time series data matrix Extract each of the clusters respectively The local median of the samples of the operation classification unit within the whole at the corresponding time point, and the global median of the samples of the overall operation classification unit at the corresponding time point, are used to generate the inter-cluster sum of squares for each time dimension. In response to the generated inter-cluster sum of squares corresponding to each time dimension and the obtained regularization parameter Generate the updated non-negative time dimension weight vector. ; In response to the updated nonnegative time dimension weight vector and the clustering partition In each of the aforementioned clusters The center point is then determined by re-identifying the actual operational classification unit sample. To update the set of center points .

[0006] Secondly, in view of the shortcomings of the prior art, the purpose of this application is to provide a computer-readable storage medium that solves the technical problems raised in the background art.

[0007] The objective of this application can be achieved through the following technical solutions: A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements a sparsity microbiome time-series data analysis method as described in the first aspect.

[0008] The beneficial effects of this application are: This invention generates inter-cluster sums of squares by extracting the local median within each cluster and the global median in each iteration. This completely eliminates the gravitational drag effect of extremely high-abundance sequencing artifacts on data centers. Based on robust feature evaluation of the median, combined with adaptive truncation of the soft threshold operator, it effectively prevents overfitting and false negative errors of the regularization penalty mechanism. This ensures that in an extremely sparse and noisy biological space, the key time windows that truly have evolutionary and differentiation significance will not be misjudged as redundant features and removed. Attached Figure Description

[0009] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0010] Figure 1This is a schematic diagram of the overall process of the sparse microbiome time series data analysis method provided in the embodiments of the present invention.

[0011] Figure 2 These are simulated feature time series diagrams under different noise levels provided in embodiments of the present invention.

[0012] Figure 3 This is a simulated data graph with an original noise level of 4.5 provided in an embodiment of the present invention.

[0013] Figure 4 This is a graph showing the results of the PAM clustering method provided in this embodiment of the invention on simulated data.

[0014] Figure 5 This is a graph showing the results of the K-means clustering method provided in this embodiment of the invention on simulated data.

[0015] Figure 6 This is a graph showing the results of the KShape clustering method provided in this embodiment of the invention on simulated data.

[0016] Figure 7 This is a graph showing the results of the HClust clustering method provided in this embodiment of the invention on simulated data.

[0017] Figure 8 This is a graph showing the results of the LASSO-PAM++ clustering method provided in this embodiment of the invention on simulated data.

[0018] Figure 9 This is a comparison chart of the adjusted RAND Index (ARI) for various clustering methods provided in the embodiments of the present invention.

[0019] Figure 10 This is a comparison of silhouette coefficients in the empirical analysis of the Moving Pictures dataset provided in this embodiment of the invention. Figure 1 .

[0020] Figure 11 This is a comparison of silhouette coefficients in the empirical analysis of the Moving Pictures dataset provided in this embodiment of the invention. Figure 2 . Detailed Implementation

[0021] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0022] Cluster analysis of microbiome time series data is crucial for revealing the evolution of host health. However, in longitudinal studies, data often contain numerous noisy time points introduced by sequencing errors or random fluctuations in biological processes, and the data itself exhibits sparsity and heavy-tailed distribution characteristics. Existing sparse clustering methods rely on arithmetic means to calculate inter-cluster differences, making them highly sensitive to extreme outliers and leading to inaccurate feature selection. This application introduces a robust inter-cluster sum of squares calculation based on the median, combined with soft threshold regularization and an improved centroid initialization strategy, to provide an analytical method that can resist outlier interference, accurately screen key temporal features, and achieve stable clustering.

[0023] Figure 1 A flowchart illustrating a method for analyzing sparse microbiome time-series data is shown. Figure 1 As shown, embodiments of the present invention provide a method for analyzing sparse microbiome time series data, the method comprising: Obtaining the time-series data matrix of the microbiome The time series data matrix Include Samples of each operational classification unit in The values ​​at each aligned time point Characterizing the time series data matrix The first in Samples of each operational classification unit; Get the set number of clusters and regularization parameters ; In response to the acquired time series data matrix and the number of clusters From the time series data matrix Extract initialization from Each center point set of center points And based on the time series data matrix The time dimension attribute initialization includes Each weight value Non-negative time dimension weight vector ; In response to the acquired time series data matrix The initial set of center points and the non-negative time dimension weight vector The algorithm iteratively solves the problem by minimizing the following objective function, determining whether the parameter changes satisfy a preset convergence condition. If the result indicates that the preset convergence condition is met, the current clustering partition is output. The set of center points and the non-negative time dimension weight vector As the optimal combination of solutions to the objective function; The objective function is: in, The characterization includes Cluster Clustering partitioning, Samples characterizing the operational classification unit To its cluster center point The weighted squared distance; The iterative solution includes: Obtain the samples of the operation classification unit To each of the aforementioned center points The weighted distance is used to classify the samples of each of the operation classification units. Uniquely assigned to the nearest center point The corresponding cluster In the process, the clustering partition is generated. ; In response to the generated cluster partition From the time series data matrix Extract each of the clusters respectively The local median of the samples of the operation classification unit within the whole at the corresponding time point, and the global median of the samples of the overall operation classification unit at the corresponding time point, are used to generate the inter-cluster sum of squares for each time dimension. In response to the generated inter-cluster sum of squares corresponding to each time dimension and the obtained regularization parameter Generate the updated non-negative time dimension weight vector. ; In response to the updated nonnegative time dimension weight vector and the clustering partition In each of the aforementioned clusters The center point is then determined by re-identifying the actual operational classification unit sample. To update the set of center points .

[0024] Specifically, we obtain a time-series data matrix X, which includes a relative abundance matrix of individual microbial sequencing data and time-node information. And the relative abundance matrix Linear mixture model spline fitting is performed to generate the time series data matrix X. This process effectively reduces high-frequency sequencing noise and inter-individual variability in the original data through model fitting, providing smoother and more representative input data for subsequent cluster analysis. As part of data preprocessing, the original sequencing data can be manipulated before fitting. For example, the original operational taxon abundance matrix and corresponding metadata table can be extracted from the original sequencing data. The original operational taxon abundance matrix can be filtered to remove operational taxonomies with a total relative abundance below a preset threshold (e.g., 0.01%). The filtered absolute counts can be converted into relative abundance to eliminate sequencing depth differences, and the sampling time in the metadata table can be standardized and mapped to a unified interval. Subsequently, the preprocessed relative abundance matrix and standardized time vector are used as input to perform linear mixture model spline fitting. During the fitting process, for each operational taxonomy unit, its relative abundance is used as the dependent variable and standardized time is used as the independent variable to construct the model. The common time trend at the population level is modeled as a fixed effect and fitted using a spline function. At the same time, the baseline differences and evolution rate differences between individuals are separated as random effects. After iteratively solving the model, only the smoothed spline curves calculated from the fixed effects parameters are extracted to eliminate individual differences and high-frequency noise, outputting the continuous fitted profile data of the operational classification unit. Finally, the fitted profile data of all operational classification units are standardized (e.g., z-score standardization) and arranged in a structured manner to generate the required time series data matrix X.

[0025] Specifically, in this application, before fitting the raw individual-level microbial time series data to LMMS, the raw individual-level microbial time series data can be preprocessed as follows: First, the abundance matrix of raw operational taxa is extracted from the acquired microbiome sequencing data. and the corresponding metadata table ; wherein, the original operational classification unit abundance matrix The elements are the absolute counts of each operational taxonomic unit in each sequencing sample, and the metadata table is... Includes a unique sample identifier, gender, and sampling time for each sequencing sample. .

[0026] For ease of understanding, this application uses data from a longitudinal microbiome study published in 2011 by Caporaso et al. Using high-throughput sequencing technology, the microbial communities of two healthy subjects (one male and one female) were sampled frequently over more than 15 months, creating the classic "Moving Pictures" dataset. This dataset, with its high sampling frequency, long time span, and rigorous data quality, has become one of the benchmark datasets for the development and validation of microbiome time-series analysis methods. This application extracted 16S rRNA gene sequencing data from gut microbial samples from this dataset as microbiome sequencing data.

[0027] Table 1 shows examples of operational taxonomic unit counts for some samples of 16S rRNA gene sequencing data from gut microbiota. In these samples, the operational taxonomic unit count is 0, which reflects the sparsity of the microbial data.

[0028] Table 1 shows examples of partial original operation classification unit count data. Specifically, the original absolute count matrix containing the specific sequence identifier (e.g., 550.L2S197.s.2.sequence) is extracted as the original operational classification unit abundance matrix. The rows of this matrix represent different operational taxonomic units, the columns represent different sequencing samples, and the elements in the matrix (e.g., 0 or 1) represent the absolute number of reads obtained from sequencing. Simultaneously, 550.L2S197.s.2.sequence extracts the corresponding sample metadata table as the metadata table. The table specifically records the gender (e.g., the host gender is set to "female" or "male") corresponding to the unique identifier of the aforementioned sample, as well as the sampling time. (For example, recording the actual number of days since the start of the experiment, such as "140 days" or "177 days"); for example, the metadata table of partial sample metadata of the 16S rRNA gene sequencing data of the gut microbiome sample exemplified in this application. As shown in Table 2.

[0029] Table 2. Metadata Tables of Some Sample Metadata Examples Secondly, the abundance matrix of the original operational classification units Filtering is performed to remove operational classification units whose total relative abundance in all samples is lower than a preset abundance threshold, generating a filtered abundance matrix. This aims to reduce random noise interference caused by extremely low-frequency sequencing errors at the source of data processing.

[0030] For the purpose of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbial sample is used as microbiome sequencing data: the abundance matrix of the original operational taxonomic units... Statistical analysis was performed to calculate the total relative abundance of each operational taxonomic unit in all intestinal samples, and the preset abundance threshold was set to... and the total relative abundance is less than or equal to The entire row of data containing the operational classification unit is removed, thereby accurately filtering out fragmented reads that are not statistically significant due to extremely shallow sequencing depth or occasional errors, and generating the filtered abundance matrix. Regarding the example in this application, the preset abundance threshold is set to In practice, other values ​​can be set according to the actual accuracy requirements.

[0031] Next, the total sequencing depth corresponding to each sequencing sample is obtained, and the filtered abundance matrix is ​​then... The absolute count values ​​in the matrix are divided by the total sequencing depth of the corresponding sample to convert them into a relative abundance matrix. This eliminates systematic baseline bias caused by uneven sequencing data between different samples, ensuring that the data between samples are objectively comparable.

[0032] For the purpose of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as the microbiome sequencing data: the filtered abundance matrix is ​​read. The total number of sequencing reads for each intestinal sample column is defined as the "total sequencing depth" of that sample. Then, the absolute count of each retained operational taxonomic unit under that sample column is divided by the total sequencing depth of that sample to obtain a standardized ratio. This division operation is performed on all sample columns to completely convert the absolute count form into a percentage form, generating the relative abundance matrix. .

[0033] At the same time, from the metadata table Extracting the discrete sampling time And by linear scaling the sampling time Mapped to a pre-defined unified time interval, a standardized time vector is generated. .

[0034] For ease of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota samples are used as microbiome sequencing data: given that the actual sampling days of human gut microbiota samples exhibit an irregular discrete distribution (e.g., day 70, day 104, day 140, etc.), from the aforementioned metadata table Extract these discrete time values ​​as the sampling time. Next, linear mapping algorithms such as range normalization are applied to map and compress the actual sampling days spanning several months to a preset closed interval (e.g., ...). Within this numerical range, generate a continuous and standardized time vector that eliminates differences in physical time scale. .

[0035] Finally, the relative abundance matrix generated by the above preprocessing is... and the corresponding standardized time vector Specifically, as input data for the subsequent LMMS fitting: For the relative abundance matrix Each operational taxonomic unit (OTU) in the dataset uses its relative abundance in different samples as the response variable. ,in Indicates an individual (sample or subject). This represents the observation time point corresponding to the individual; expressed as a standardized time vector. A linear mixed-effects spline model is constructed using these as independent variables. Specifically, the original time variables are first extracted from the metadata and then mapped to a unified time interval through linear normalization: A mapping from continuous-time to discrete-time grid {1,…,60} is implemented to eliminate the influence of irregular sampling. Simultaneously, sample identifiers are introduced into the model as individual variables to characterize the repeated measurement structure. Subsequently, the time effect is expanded based on P-spline basis functions, and the nonlinear time trend at the population level is represented by the following linear mixed-effects model: in For the fixed effects portion, and Each represents an individual The random intercept and random slope satisfy: And with the error term They are independent of each other.

[0036] During model fitting, a smoothing parameter is introduced through penalized likelihood estimation. For fixed effects coefficient Apply a difference penalty and estimate using the restricted maximum likelihood method: in Indicates the first Each spline basis function (usually a B-spline basis). This is the corresponding fixed effects coefficient. The number of basis functions (determined by the number of nodes or automatically selected by the algorithm); The difference order (usually taken as) This function controls the smoothness of the curve. In practice, it can be called directly. lmmSpline(data,time,sampleID,basis="p-spline",timePredict=1:60) The "sample × OTU" matrix is ​​modeled column-by-column; after model convergence, only the smoothed spline predictions corresponding to the fixed effects are extracted (i.e., predSpline output, corresponding to the code modeled.data<-t(slot(lmms.output,'predSpline'))), thereby removing the influence of individual random perturbations and sequencing noise; finally, the predictions of all OTUs on a unified time grid are combined to form a smoothed expression matrix. The data is then further converted into long-format data (melt(modelled.data_f) in the code) for subsequent time series clustering analysis and dynamic pattern recognition.

[0037] For ease of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota samples are used as microbiome sequencing data: for a specific bacterium in the gut, spline interpolation technology is used to complete its numerical distribution over the unsampled days; at the same time, since there are physical individual differences in the basic gut microbiota of the two healthy subjects (one male and one female), the system removes this baseline difference caused by gender and physical condition as a "random effect" at the data level, and only extracts the fitted contour data. To characterize the smooth, common trajectory curves of this type of bacteria's evolution over time in the human gut.

[0038] Finally, the fitted contour data Perform a traversal, remove columns with all zero values, and generate non-zero profile data. Subsequently, the non-zero contour data Perform data standardization processing to scale the abundance curves of each operational taxonomic unit to a uniform dimensionless scale, generating standardized data. Next, the standardized data Convert to long format data .

[0039] After the above standardization process, the retained Each operational classification unit in The values ​​at each aligned time point are arranged in a structured manner to generate a time series data matrix. .

[0040] For ease of understanding, the 16S rRNA gene sequencing data of the aforementioned gut microbiome samples are used as microbiome sequencing data. Thousands of smoothed gut microbiome evolution curves are arranged in rows, with each column strictly corresponding to an aligned, standardized time point. The absolute abundance values ​​of all curves are uniformly scaled to the same evaluation scale (e.g., mean 0, variance 1), completely eliminating the order-of-magnitude barrier between high-abundance and low-abundance bacteria. This ultimately generates and solidifies a microbiome containing... OK The structured numerical matrix of columns is the time series data matrix. Thus, the time series data matrix is ​​complete. Enter the information in step 2.

[0041] LMMS was used to fit the original microbial time series data at the individual level to generate fitted time series profile data. The algorithm is implemented as follows: (2) Step 2: PAM++ initializes cluster centers The PAM++ initialization strategy, based on the K-means++ idea, is adopted to gradually select from the dataset. Initialize medoids. A distance-sensitive probability sampling mechanism is used to distribute the initial medoids as widely as possible in the weighted space, thus avoiding the instability caused by random initialization. The initialization formula is as follows: in This is the set of currently selected center points. This is the weight matrix.

[0042] The specific implementation is as follows: S201. First, obtain the time series data matrix. The set number of clusters and non-negative time dimension weight vector Based on the nonnegative time dimension weight vector Constructing the weight matrix Subsequently, random numbers from the time series data matrix were obtained using a uniformly distributed random number generator. Included The first sample is randomly selected from the samples of each operational classification unit as the first medoid. And store it in the currently selected center point set. middle; For the purpose of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as microbiome sequencing data: assuming the system is pre-defined to divide gut microbiota into... There are several evolutionary groups. First, an initial weight representing the importance of each sampling day is extracted using a uniform initialization method, meaning each sampling day has an equal weight. A diagonal weight matrix is ​​then constructed with these weights as its diagonal elements. Next, from the set of operational taxonomic units containing thousands of gut microbes, a microbe is selected by random sampling instructions (e.g., the evolution curve corresponding to "550.L2S197.s.2.sequence") and stored in the set as the center point of the first taxonomic group.

[0043] S202. Let the set of currently selected center points be... The weighted minimum distance is defined as: in It is a diagonal matrix weighted by the feature dimensions. The distance-weighted probability sampling formula is: For each operational classification unit that has not yet been selected, the sample In response to the weight matrix and the current set of selected center points Calculate the samples respectively To the set Add to all existing center points The minimum value of the squared weighted Euclidean distance is taken as the weighted minimum distance. The purpose is to introduce a weight matrix of the time dimension in advance when calculating multidimensional spatial distance, so as to suppress the high noise time points in the initial stage and highlight the feature differentiation of key time intervals.

[0044] For ease of understanding, the 16S rRNA gene sequencing data of the aforementioned gut microbial samples are used as microbiome sequencing data: For each remaining microorganism in the gut sample that has not yet been selected as a centroid, the spatial difference between its evolution curve and the curve of the microorganism currently selected as a centroid is calculated using a weighted distance formula. If there are multiple centroids in the current set, the distance between the microorganism and the nearest centroid is extracted as its weighted minimum distance. .

[0045] S203, Obtain the weighted minimum distance of all candidate operation classification units. Then, probability normalization is performed. Specifically, the weighted minimum distance sum of all samples of the operational classification unit within the entire set is calculated. And according to the formula Convert the absolute distance values ​​of each sample into a range of The normalized probability distribution within the medoid is used. At this point, the further the operation classification unit is spatially from the existing center point, the higher its selection probability. Subsequently, weighted random sampling is performed based on this normalized probability distribution to select the next operation classification unit as the new medoid. and add it to the set of currently selected center points. middle.

[0046] For ease of understanding, as an example, the 16S rRNA gene sequencing data of the aforementioned gut microbiome sample is used as the microbiome sequencing data: the weighted minimum distances of all remaining microorganisms calculated above are summed, and then the proportion of each microorganism's weighted minimum distance to the sum is calculated as its probability value for selection. This means that the new species that differs more significantly from the selected gut microorganisms in its dynamic evolutionary pattern has a higher probability of being selected as the next new taxa center. Random sampling instructions are strictly executed according to the set probability distribution to extract the next microorganism to be added to the center set.

[0047] S204. Repeat S201-S203 until the number of elements in the currently selected center point set reaches the set number of clusters. At this point, the loop terminates, and the corresponding initialized set of center points is output for the latest set of center points. .

[0048] In this way, the representative points in the initialized set of centroids are maximized in the data space, thus providing a highly robust topological starting point for subsequent clustering and allocation.

[0049] For ease of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbial samples are used as microbiome sequencing data: the above process is repeated continuously until six gut microbial evolution curves with significant evolutionary differences are successfully selected. These six curves, selected through probability sampling, constitute the initialized set of center points. This effectively avoids cluster collapse caused by multiple initial center points representing the same type of similar microorganisms due to purely random sampling.

[0050] In response to the acquired time-series data matrix X, the initialized centroid set M, and the non-negative time-dimensional weight vector w, an iterative solution is performed with the constraint of minimizing the following objective function J, wherein the objective function J is: in, The characterization includes Cluster Clustering partitioning, Samples characterizing the operational classification unit To its cluster center point The objective function is defined by the weighted squared distances; λ is the regularization parameter. The first term measures the internal compactness of the cluster partition, i.e., the sum of the weighted distances from all samples to their respective cluster centers. The second term is an L1 norm penalty term for all weights, the strength of which is controlled by the regularization parameter λ. Its function is to drive the time dimension weights that contribute little to the clustering towards zero, thus achieving feature selection. The constraints ensure that the weight vector is a valid probability distribution. The iterative solution process is to find the optimal combination of the cluster partition C, the set of center points M, and the non-negative time dimension weight vector w that minimizes the objective function J.

[0051] In this embodiment, the iterative solution process includes steps 3-7, wherein: (3) Step 3: Sample allocation Given the current weight vector and the set of center points Under the condition of weighted Euclidean distance, Each sample is assigned to the cluster corresponding to the nearest medoid, resulting in a new clustering partition. The weighted distance calculation formula is: The specific implementation process is as follows: S301. First, obtain the current non-negative time dimension weight vector. and the current set of center points For the time series data matrix Samples of each operational classification unit in the process Calculate the sample respectively With the set of center points Each center point spatial distance , The formula is as follows: in, and They represent the samples to be assigned. With the center point In the Standardized abundance values ​​at each time point Indicates the first Each time point is assigned a non-negative weight. The purpose is to introduce weights when calculating the overall evolutionary morphological differences between samples. By linearly scaling the squared error at each time point, the dominant role of key time intervals in allocation decisions is amplified in terms of spatial metrics, while suppressing interference from irrelevant or noisy time intervals.

[0052] For ease of understanding, this application is illustrated by using the 16S rRNA gene sequencing data of the aforementioned gut microbiome samples as microbiome sequencing data: For a specific microorganism in the gut sample that has not yet been classified (i.e., the operational taxonomic unit to be assigned), its evolution curve is compared with the curves of the six currently selected representative microorganisms (i.e., six medoids) for similarity. When comparing the abundance difference between two time points (e.g., day 140), the absolute difference is not used directly; instead, the square of the difference is multiplied by the specific weight currently assigned to the "day 140" time point by the system. If the data from a certain number of days contains high sequencing noise, its corresponding weight... The difference in abundance over the past few days will be so small that it will have almost no impact on the final total distance, thus ensuring the anti-interference ability of the similarity comparison.

[0053] S302, After calculating the sample To all After calculating the weighted Euclidean distance of each center point, the minimum distance value is extracted, and this sample is... It is uniquely assigned to the cluster corresponding to the medoid that produced the minimum distance.

[0054] S303, regarding the time series data matrix All of them Each operation classification unit executes steps S301-S302 in a loop. After all samples have been traversed, the global clustering partition is generated and output. This clustering partition The cluster affiliation of each operational classification unit under the current parameter conditions is clearly recorded and will be directly passed to subsequent steps as the basic data structure to evaluate the internal density of each cluster and the differences between clusters.

[0055] For ease of understanding, the 16S rRNA gene sequencing data of the aforementioned gut microbial samples are used as microbiome sequencing data. After calculating the weighted distance between the specific microorganism and six representative microorganisms, the representative with the smallest distance (i.e., the most similar dynamic evolutionary pattern) is identified, and this specific microorganism is directly labeled with the cluster to which that representative belongs. This calculation and labeling process is performed on thousands of gut microorganisms one by one, ultimately clearly dividing this vast and chaotic microbial community into six non-overlapping evolutionary groups, generating the current clustering set. .

[0056] (4) Step 4: BCSS calculation Under the current clustering partition, calculate the inter-cluster sum of squares (BCSS) for each time dimension, and use the median to enhance robustness: in and They represent the first Clusters and the overall sample in the first The median at each time point.

[0057] After obtaining the current complete clustering partition, it is necessary to independently evaluate the distributional dissimilarity between clusters for each time dimension, i.e., calculate the Between-Cluster Sum of Squares (BCSS). The specific implementation process is as follows: S401, Clustering partitioning in response to S303 and the input time series data matrix For time series data matrices The At the first time point, the standardized abundance values ​​of all operational taxonomic units at that time point are first sorted, and the overall sample is extracted at the second time point. median at each time point Simultaneously, for each individual cluster Count the number of operational classification units contained in the cluster. And extract all samples within that cluster in the first... median at each time point The aim is to abandon the traditional arithmetic mean and adopt the "median" as the statistical benchmark, which can greatly enhance the robustness of the algorithm and avoid statistical bias caused by individual extreme abnormal operation units (such as noisy species with sudden over-proliferation) at a single point in time.

[0058] For the purpose of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as microbiome sequencing data: assuming that the evaluation is currently being conducted on "day 140 (i.e., day...)". The data performance at each time point was analyzed. First, the abundance values ​​of all thousands of gut microbiota at day 140 were sorted, and the middle value was taken as the "global median". Then, the data for each group (cluster) was analyzed separately. The study analyzed hundreds of microorganisms and identified the "local median" of their abundance on day 140, recording the total number of microorganisms in that group. This median-based mechanism ensures that even if a specific gut bacterium experiences an "overflowing" abundance on day 140 due to sequencing anomalies or sudden diarrhea in the host, it will not skew the statistical representativeness of the entire group.

[0059] S402, Calculate each cluster Local median Compared with the global median Square the difference between them, and multiply the squared value by the number of samples in the corresponding cluster. and all The product results of each cluster are summed to obtain the result. , The formula is as follows: S403, for all Execute at each time point. The calculation generates a vector set of inter-cluster sums of squares containing the variability assessment values ​​for each time dimension. The numerical value intuitively quantifies the first The degree of dispersion of the evolution patterns among different clusters at each time point; the larger the value, the higher the contribution of that time point to distinguishing different cluster structures.

[0060] For ease of understanding, this application is illustrated by using the 16S rRNA gene sequencing data of the aforementioned gut microbiome sample as microbiome sequencing data: For the aforementioned "day 140," the difference between the median of the first group and the global median is calculated, squared, and multiplied by the number of microorganisms in that group. The values ​​for the other five groups are calculated using the same method, and all six results are summed to obtain the "BCSS value for day 140." If this value is very large, it indicates that the six gut microbiome groups exhibit distinctly different behaviors on day 140, meaning this day is a "critical identification period" for distinguishing them; conversely, if the value is small, it indicates that all microorganisms behave similarly on this day, lacking differentiation. This "criticality score" for each day is calculated and directly fed into the next step, the LASSO regularization module, for weight cleaning.

[0061] (5) Step 5: Weight update based on LASSO The time dimension weights are updated based on the BCSS values, and a soft threshold operator is introduced to achieve LASSO sparsity: This step can automatically filter key time intervals that significantly contribute to the cluster structure.

[0062] The specific implementation process is as follows: S501, Obtain the set LASSO regularization parameters And the inter-cluster sum of squares for each time dimension in S403 Calculate the inter-cluster sum of squares With the LASSO regularization parameter The difference is calculated, and a conditional judgment instruction is executed. If the difference is greater than zero, the difference is retained as the valid feature output value at that time point; if the difference is less than or equal to zero, the valid feature output value at that time point is forcibly truncated to zero. The purpose is to apply a uniform numerical contraction to the feature evaluation values ​​of each time dimension by introducing a penalty term, so that redundant time points or high-noise time points that contribute little to clustering and have insignificant differences between clusters are completely eliminated, thereby suppressing the interference of high-frequency noise time points on the clustering results.

[0063] For ease of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as microbiome sequencing data: The set regularization parameters are obtained. Subsequently, the inter-cluster sum of squares score corresponding to each previously calculated sampling day (e.g., day 70, day 140) is extracted sequentially. If the score for a certain sampling day is less than or equal to the regularization parameter... If the differences in microbial evolution among different groups at this time point do not reach a significant level, a soft threshold truncation mechanism is directly triggered, forcibly truncating the feature score corresponding to the sampling day to zero; if the score of the sampling day is greater than the regularization parameter... If the value exceeds the regularization parameter, then the value is calculated and retained.

[0064] S502. After completing the S501 operation, sum the non-zero valid feature output values ​​retained throughout the entire time period to obtain the total sum of global valid features. Subsequently, the effective feature output values ​​at each time point after soft thresholding are divided by the sum of the global effective features to generate the updated time dimension weights for each time point. , The formula for generating the formula is as follows: S503, all of them The new weight combinations constitute the updated non-negative time dimension weight vector. After this operation, the time points with zero values ​​after processing are given absolute sparse weights of zero, while the key time points that contribute significantly to the clustering structure are reassigned non-negative weights greater than zero. Moreover, the sum of the weights of the entire time dimension is strictly equal to 1, thus realizing dynamic dimensionality reduction and key dimension purification of high-dimensional time series features.

[0065] For ease of understanding, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as microbiome sequencing data: The scores of all sampling days with values ​​greater than zero after truncation are summed to calculate the total sum of global effective features. Then, the effective score retained for each sampling day is divided by this total sum of global effective features to obtain the updated weight ratio for each sampling day (i.e., the updated time dimension weight). The sampling days whose scores were truncated to zero had strictly zero new weights, meaning the system would not include sequencing data from these specific sampling days in subsequent distance metric and clustering iterations; while the sampling days that retained valid scores were proportionally allocated a total of The weighting of the data is calculated. This operation enables the automatic filtering out of invalid sequencing sampling points in the process of gut microbiome evolution, allowing the system's clustering computing power to be precisely focused on key sampling nodes with high discriminative power.

[0066] (6) Step 6: Medoid update With the updated weights and clustering partitions, a new medoid is selected for each cluster to minimize the sum of weighted squared distances within the cluster: After obtaining the updated feature weights, the center points of each cluster are re-evaluated and iteratively updated based on the purified high-dimensional feature space. The specific implementation process is as follows: S601, Obtain the updated non-negative time dimension weight vector after S503. and the current clustering division For each independent cluster in the clustering partitioning Iterate through all samples of all operational classification units contained within the cluster. For this cluster... Any sample to be evaluated within Based on the updated non-negative time dimension weight vector Calculate the values ​​of the sample to be evaluated respectively. With all other samples in the same cluster The weighted Euclidean distance between These squared distance values ​​are then summed to generate the sample to be evaluated. The sum of weighted squared distances within the cluster The aim is to directly apply updated weights that have been truncated by a soft threshold during the measurement process, so that time points that are judged to be noisy or redundant no longer contribute any numerical value in the distance measurement, thereby accurately measuring the local topological correlation between samples in a highly pure feature subspace.

[0067] For ease of understanding, this application is illustrated by using the 16S rRNA gene sequencing data of the aforementioned gut microbial samples as microbiome sequencing data: Among hundreds of gut microbes classified into a specific evolutionary group, one microbe is selected sequentially as the sample to be evaluated. Next, the sum of distances between the sample to be evaluated and all other microbes within that group is calculated. During the calculation, the system strictly applies the sampling day weight updated in the previous stage (i.e., sampling days whose values ​​are truncated to zero are not included in the distance accumulation, and the number of days retaining valid feature output values ​​are used for distance calculation according to their new weight ratio). The above calculation operation is performed on each microbe within the group to obtain the sum of its intra-cluster weighted squared distances.

[0068] S602, Obtaining a single cluster After summing the intra-cluster weighted squared distances of samples from all operational classification units, the sample from the operational classification unit that produces the minimum intra-cluster weighted squared distance is extracted and established as the updated medoid for that cluster, denoted as . . The formula for generating the formula is as follows: S603, for all Each cluster executes S603 in a loop, thereby combining to generate a globally updated set of center points. The aim is to require that cluster centers be actual observed samples that exist in the dataset, rather than virtual coordinate points calculated through numerical averaging; thereby directly cutting off the gravitational drag effect of extreme outliers of a single individual on the cluster center, and greatly improving the robustness of the algorithm to sudden abnormal abundance fluctuations in microbial communities.

[0069] For ease of understanding, as an example, the 16S rRNA gene sequencing data of the aforementioned gut microbial samples are used as microbiome sequencing data. By comparing the distance assessment results of all microorganisms in this specific evolutionary group, the gut microorganism with the smallest sum of weighted squared distances within the cluster (i.e., the one most closely positioned at the topological center of the group on the evolutionary curve) is identified and established as the new representative center point for this group. This mechanism of strictly selecting real microbial sequences as center points ensures that even if several bacteria with extremely abnormal abundance due to host dietary mutations are mixed into the group, the center point of the group remains a true biological curve that conforms to the common evolutionary laws of the group, and is not stretched into a virtual average curve that does not exist biologically by outliers.

[0070] (7) Step 7: Convergence determination If the changes in clustering partitioning, medoids, or weight vectors in two consecutive iterations are less than a preset threshold or reaching the maximum number of iterations. If the algorithm succeeds, it is considered converged; otherwise, return to step 3 to continue iterating. The convergence condition is: in This indicates the amount of change in medoids. This indicates the change in cluster partitioning.

[0071] The specific implementation process is as follows: S701. Obtain the non-negative time dimension weight vector after this round of iteration update. Updated set of center points And the updated clustering partitions, and retrieve the non-negative time dimension weight vector saved from the previous iteration. The previous set of centroids and the previous clustering partitions. Three core state change indicators are calculated respectively: First, calculate the change in cluster partitioning. That is, to count the number or proportion of operational classification units whose cluster affiliation has changed compared to the previous round; Secondly, calculate the change in medoids. That is, to assess the spatial changes or replacements of the set of center points established in this round compared to the set of center points in the previous round; Third, calculate the change in the non-negative time dimension weight vector. That is, calculate the absolute value of the difference between the corresponding elements of the two weight vectors at each time point, and extract the maximum absolute value (i.e., the infinite norm distance).

[0072] The purpose is to analyze the topology (including...) and The dynamic fluctuations of the algorithm are fully quantified by three independent evaluation dimensions: the clustering model's stability evaluation, the feature space (including the change in weight vectors), and the feature space (including the change in weight vectors), ensuring that the stability evaluation of the clustering model is rigorous and comprehensive.

[0073] For ease of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbial samples are used as microbiome sequencing data: the newly generated classification results are compared with the previous round of classification results to count how many gut microbes changed their taxonomic labels in this round (i.e., to calculate...). ); compare whether species replacement has occurred at the representative centers of the current six evolutionary groups (i.e., calculate) Simultaneously, compare the difference between the old and new weights for each sampling day (e.g., day 140), and find the data change on the day with the largest difference (i.e., calculate the change in weight). ).

[0074] S702. Extract the maximum value among the three change indicators mentioned above, that is... Then, the maximum value is compared with a preset threshold. Perform numerical comparisons; simultaneously, check whether the current number of iterations has reached the preset maximum number of iterations. ; If the judgment condition is met Or the current number of iterations has reached the maximum number of iterations. If the above convergence condition is not met, the loop iteration instruction is terminated; otherwise, if the above convergence condition is not met, the updated parameters of this round replace the corresponding parameters of the previous round, and return to the sample allocation operation with the updated non-negative time dimension weight vector and the updated set of centroids (i.e., jump back to step 3) to start the next round of clustering iteration optimization.

[0075] This effectively avoids the waste of computing power that leads to an infinite loop, and optimizes the scheduling of machine computing resources while ensuring high-dimensional clustering accuracy.

[0076] For ease of understanding this application, for example, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as the microbiome sequencing data: find the largest value among the three statistical variables mentioned above, and compare it with the minimum tolerance value set by the system (e.g., 0.001 as a preset threshold). The system compares the data. If the largest change is less than 0.001, it means that the classification structure of these tens of thousands of gut microbes has become fixed and no longer makes substantial adjustments, or the system has repeated the calculation to the set upper limit (e.g., 100 times), and will immediately stop the loop. Conversely, if the change is still large, it means that the classification of microorganisms is still in a period of drastic adjustment, and will overwrite the old data, take the latest calculated sampling weights for each day and the latest 6 group centroids, and return to step 3 of "comparing microbial distances and assigning groups" to start a new round of classification iteration.

[0077] (8) Step 8: Output clustering results The algorithm ultimately outputs the cluster labels of the classification units. medoids corresponding to each cluster and non-negative time dimension weight vector This enables the grouping and identification of microbial time dynamic patterns and the synchronous selection of key time intervals.

[0078] The specific implementation process is as follows: S801, Extract the final clustering partition from the output of S70. The final set of centroids corresponding to each cluster (i.e., medoids) and the final non-negative time dimension weight vector .

[0079] In practice, these three types of data can be correlated and mapped to generate a structured result list or data report containing operational classification unit identifiers, cluster labels, centroid identifiers, and time point weight ratios, which can then be output to external storage media or interactive terminals. The purpose is to output clustering partitions. It provides accurate topological attribution of samples in high-dimensional space; Output center point set It provides the most representative real evolutionary trajectories for each cluster; Output non-negative time dimension weight vector This provides a quantitative evaluation of the importance of the time dimension. Dimensions with strictly zero values ​​correspond to filtered noise or redundant time points, while dimensions with values ​​greater than zero correspond to key time points that significantly contribute to distinguishing cluster structures.

[0080] For the purpose of understanding this application, the 16S rRNA gene sequencing data of the aforementioned gut microbiota sample is used as an example, for instance: After determining that the classification process is complete and no longer subject to change, the final results are summarized into a comprehensive data report for output. In this report: 1. It records in detail the specific evolutionary group labels (e.g., labeled "Group 1" or "Group 6") to which thousands of gut microbes (operational taxonomic units) were ultimately classified; 2. The most representative real gut microbial sequences (i.e., the final center points) corresponding to each of these 6 evolutionary groups were clearly identified; 3. At the same time, a complete importance assessment table of sampling days was output, clearly showing which sampling days had the highest weight and which sampling days had their weight truncated to zero and were therefore not of analytical value.

[0081] That is, this application relies on the output clustering partitioning With the set of center points This enables highly accurate grouping and identification of massive amounts of microbial temporal dynamic patterns containing individual differences, allowing species with similar evolutionary patterns to be reliably classified; simultaneously, relying on the output non-negative time dimension weight vector... It enables adaptive selection and extraction of key time intervals, thereby allowing the computational model to automatically lock the core time window that causes differences in community structure in extremely sparse and noisy longitudinal microbiome data without human intervention.

[0082] Verification 1: To verify the effectiveness of the LASSO-PAM++ algorithm, the following simulation experiments were conducted: Introduction to Simulation Data In the simulation experiment analysis section, this application mainly refers to the simulation framework of Kodikara et al. to generate representative microbiome time series data, and designs experimental comparisons for different dynamic modes, noise levels and sample structures.

[0083] First, this application generated simulated time series The time point is from arrive And distribute these sequences evenly to Each group, each class A sequence. (Note: The original text contains several characters and symbols, which are difficult to translate accurately without further context.) Indicates the first A set of simulated sequences for each class group.

[0084] Subsequently, for each simulation sequence This application, at each time point Add independent and identically distributed Gaussian noise (Representing variation between individuals) to simulate from Observational sequences of different individuals: in Indicates the number of time points. This indicates the noise level. This simulation strategy is consistent with previous studies on longitudinal microbiome simulations, which use Gaussian noise to simulate periodic individual differences to reflect variability among real observations.

[0085] For each noise level (This application sets three noise levels) (Used to examine the impact of inter-individual variation on clustering performance), this application generates for each simulated sequence. simulated sequence Thus, a containing A simulated sample set of time series. To improve statistical robustness, samples are repeatedly generated at each noise level. Independent datasets ( This simulation design references existing longitudinal microbiome clustering evaluation frameworks, which use noise level as a key parameter to quantify the clustering algorithm's ability to identify real dynamic patterns.

[0086] For each simulation dataset This application first applies a linear mixture model spline fitting (LMMS) to the simulated sequences for smoothing, in order to reduce noise interference and estimate the continuous trajectory of each sequence. This represents the fitted spline function value, i.e.: in This represents the fitting function based on linear mixture model splines, which can capture inter-individual differences and time dependence by embedding random effects. The trajectory matrix obtained after LMMS fitting. It is used as input to clustering algorithms.

[0087] Evaluation indicators To comprehensively evaluate the performance of the proposed algorithm, this application selected four commonly used clustering evaluation metrics: the Rand Index (RI), the Adjusted Rand Index (ARI), the Jaccard Index (JI), and the Normalized Mutual Information (NMI). These metrics measure clustering effectiveness from different dimensions, providing a quantitative basis for comparing algorithms.

[0088] (1) RAND Index The RAND Corporation (RI) measures the performance of an algorithm by comparing the clustering results with the true class. The RI value ranges from 0 to 1, with values ​​closer to 1 indicating greater consistency between the clustering results and the true class. The formula is as follows: in, This represents the number of sample pairs that belong to the same class in both the clustering results and the true class. This represents the number of sample pairs that belong to different classes in both the clustering results and the true class. This represents the number of sample pairs that belong to the same cluster in the clustering results but belong to different clusters in the actual clustering. This represents the number of sample pairs that belong to different classes in the clustering results but belong to the same class in the actual categories.

[0089] (2) Adjust the RAND index The Adjusted Rand Index (ARI) introduces a correction term to the Rand Index, eliminating the influence of random clustering on the evaluation of results. A higher ARI value indicates better clustering results. The formula is as follows: in, This represents the original RAND index; The expected value of the RAND index is calculated under the condition of random clustering. The maximum value of the RAND index is calculated under the condition of optimal clustering.

[0090] (3) Jaccard Index The Jaccard Index (JI) measures similarity by comparing the clustering results with the proportion of identical samples in the true cluster. The JI value ranges from 0 to 1, with higher values ​​indicating better clustering. The formula is as follows: in, , , The definition is the same as in the RAND Corporation index.

[0091] (4) Normalized mutual information Normalized mutual information (NMI) measures the degree of information sharing between the clustering results and the true clusters. NMI values ​​range from 0 to 1, with higher values ​​indicating more information shared between the clustering results and the true clusters. The formula for mutual information (MI) is as follows: in, yes and Simultaneously take values and The joint probability distribution, and They are and The marginal probability distribution. The NMI formula is as follows: in, It is the clustering result and real category labels Mutual information between them; It is the clustering result Entropy; These are real category labels. The entropy.

[0092] Visualization and Comparative Analysis To comprehensively evaluate the performance of the proposed LASSO-PAM++ algorithm in time series clustering tasks, a time series dataset was generated using the simulation described above, and a systematic clustering algorithm performance evaluation experiment was conducted. This dataset is generated using four basic patterns ( Figure 3 The dataset (original Noise=4.5 data) was generated, containing 200 time series of length 50, divided into 4 categories. This application compares PAM, K-means, K-Shape, HClust, and the LASSO-PAM++ algorithm proposed in this application. Figures 4-8 The clustering results of each algorithm are visualized. Figure 9 Quantitative performance comparisons were then conducted by adjusting the RAND Index (ARI).

[0093] In the visual comparison experiment, the LASSO-PAM++ algorithm demonstrated excellent performance. For example... Figure 8 As shown, the clustering results generated by this algorithm exhibit a clear cluster distribution: the data distribution is clear, the cluster boundaries are obvious, and there is less confusion between time series data of different categories. Compared with other algorithms, LASSO-PAM++ has better separation of time series patterns and shape features in different categories, effectively avoiding the problem of inter-class overlap in common algorithms. Meanwhile, within the same category, the time series show high compactness, with data with similar patterns distributed more closely together, demonstrating good intra-class consistency. These advantages are quantitatively verified in the ARI comparison chart, further confirming the effectiveness and superiority of the LASSO-PAM++ algorithm. (Figure...) Figure 3 This is a graph of the original Noise=4.5 data; Figure 4 The results are from the PAM clustering analysis. Figure 5 The results are from the K-means clustering. Figure 6 The result is from KShape clustering. Figure 7 The results are from the HClust clustering. Figure 8 The results are from LASSO-PAM++ clustering. Figure 9 This is a comparison chart of clustering ARI.

[0094] Comparison of data metrics for clustering algorithms To evaluate the clustering performance of K-means, K-medoids, HAC, K-shape, and the proposed LASSO-PAM++ algorithm under different noise levels, this section conducts experiments on the simulated datasets RDATA15 and RDATA45. To eliminate random errors, each clustering method is tested 100 times, and the average values ​​of four evaluation metrics—RI, ARI, NMI, and JI—are calculated to systematically compare and analyze the performance of each algorithm.

[0095] Table 3. Performance comparison of different clustering algorithms on the RDATA15 dataset (noise level 1.5) Table 3 shows the performance of the five clustering algorithms under low-noise conditions (RDATA15, noise level 1.5). The LASSO-PAM++ algorithm achieves the best results in all four evaluation metrics, demonstrating its leading performance advantage.

[0096] The Adjusted Rand Index (ARI) is used to measure the consistency between clustering results and true labels. Its value range is [-1, 1], with a higher value indicating better clustering performance. On the ARI metric, the LASSO-PAM++ algorithm achieved the highest score of 0.874, a 4.92% improvement over PAM clustering (0.833).

[0097] Normalized mutual information (NMI) is used to evaluate the degree of information sharing between clustering results and true labels, with a value range of [0,1]. The LASSO-PAM++ algorithm achieved the highest NMI score of 0.870, a 5.20% improvement over PAM clustering (0.827).

[0098] The RAND Index (RI) measures the similarity between cluster assignments and true labels, with a value range of [0,1]. On the RI index, the LASSO-PAM++ algorithm achieved the highest score of 0.953; K-means (0.943) and PAM (0.938) also performed well, but LASSO-PAM++ remains the leader.

[0099] The Jaccard Index (JI) assesses the similarity and accuracy of clustering results, with a value range of [0,1]. The LASSO-PAM++ algorithm achieved the highest JI score of 0.826, a 6.44% improvement over PAM clustering (0.776).

[0100] Based on four evaluation metrics, the LASSO-PAM++ algorithm outperforms other comparative algorithms under low-noise conditions and exhibits high stability with zero standard deviation in 100 repeated experiments. This algorithm achieves feature selection and sparse representation through LASSO regularization, combined with an improved PAM++ cluster center initialization strategy, effectively enhancing the clustering accuracy and robustness of microbiome time series data.

[0101] Table 4. Performance comparison of different clustering algorithms on the RDATA45 dataset (noise level 4.5) In terms of the ARI (Average Recognition Index), as shown in Table 4 above, the LASSO-PAM++ algorithm achieved the highest score of 0.837, outperforming other comparative algorithms. Compared with K-shape clustering (ARI=0.787), LASSO-PAM++'s ARI value improved by 6.35%, demonstrating its advantage. This result indicates that even under high-noise conditions, LASSO-PAM++ can still accurately identify the true clustering structure of microbiome time series data, and its sparsity and regularization characteristics help to more accurately capture the inherent grouping patterns of the data.

[0102] On the NMI metric, the value ranges from 0 to 1, with a higher value indicating a greater degree of information sharing between the clustering results and the true category labels. The LASSO-PAM++ algorithm achieved the highest score of 0.838, outperforming other comparative algorithms. This demonstrates that LASSO-PAM++ excels in capturing the category information of microbiome time series data under high-noise environments, and its clustering strategy combined with LASSO sparsity penalty can better identify the potential structural features of the data.

[0103] On the RI metric, the value ranges from 0 to 1, with a higher value indicating a closer match between the clustering results and the true class labels. The LASSO-PAM++ algorithm achieved an excellent score of 0.939, outperforming other comparative algorithms. This result fully demonstrates the superior performance of LASSO-PAM++ in accurately assigning microbiome time-series data points to their corresponding classes under high-noise conditions. Its improved cluster center selection mechanism based on PAM++ enhances the algorithm's accuracy and stability.

[0104] The Jaccard index (JI) is used to measure the similarity and accuracy of clustering algorithm results. The JI value ranges from 0 to 1, with a higher value indicating a higher similarity between the clustering results and the true class labels. The LASSO-PAM++ algorithm achieved the highest score of 0.826, outperforming all other compared algorithms. Compared to PAM clustering (JI=0.741), it represents an improvement of 11.47%. This result demonstrates that LASSO-PAM++ can still produce grouping results highly similar to the true labels even in noisy environments, and its ability to combine LASSO feature selection helps improve the purity and consistency of the clustering results.

[0105] Comprehensive performance analysis Based on the combined results of the four clustering performance evaluation metrics, LASSO-PAM++ achieved the best performance in all metrics under both low and high noise conditions. It comprehensively outperformed traditional PAM clustering, K-means clustering, hierarchical clustering, and K-shape clustering methods in all four dimensions: ARI, NMI, RI, and JI. This indicates that LASSO-PAM++ has higher accuracy and robustness in clustering results on microbiome time series data.

[0106] The superior performance of the LASSO-PAM++ algorithm is mainly attributed to the following two innovations: First, by introducing LASSO regularization, feature selection and sparse representation are achieved, effectively removing high noise and redundant information in microbiome time series data; second, the improved PAM++ cluster center initialization strategy enhances the algorithm's convergence speed and clustering quality. This combined advantage makes LASSO-PAM++ particularly suitable for clustering tasks of microbiome time series data with high dimensionality, temporal correlation, and high noise characteristics.

[0107] Verification 2: To objectively evaluate the quality of microbiome time-series clustering results, this application selected the silhouette coefficient as an indicator of internal cluster validity. Originally proposed by Rousseeuw, this index measures cluster tightness and inter-cluster separation and has been widely applied in cluster analysis theory and application research. In microbiome research, the silhouette method is also used to select the optimal number of clusters and evaluate the effectiveness of community partitioning. In a human microbiome community structure analysis published in Genome Biology, Huttenhower et al. used the silhouette method based on 16S rRNA data to determine the optimal number of clusters for microbial communities in different body parts and used the average silhouette value to judge the rationality of the clustering results.

[0108] For each sample ,set up This is the average distance between the sample and other samples within its cluster. The silhouette coefficient of a single sample is defined as the average distance between the sample and all samples in its nearest neighbor clusters. The quality of the overall clustering results is usually assessed using all samples. The mean silhouette coefficient is used to measure the sample size. The silhouette coefficient ranges from -1 to 1, where a value closer to 1 indicates that the samples within a cluster are more closely packed and the samples between clusters are more separated; a value closer to 0 indicates that the sample is located at the cluster boundary; and a negative value means that the sample may be misclassified into the current cluster.

[0109] Empirical Results and Analysis Algorithm superiority: Quantization comparison based on profile coefficients like Figure 10 and Figure 11 In the high-resolution microbiome time-series dataset "Moving Pictures of the Human Microbiome," which includes two subjects, four body parts, and 396 time points, LASSO-PAM++ significantly outperforms traditional clustering methods in the internal clustering quality metric—the silhouette score. Compared to the suboptimal K-means (0.807) and Ward hierarchical clustering (0.799), LASSO-PAM++ achieves the highest silhouette score of 0.876, representing a substantial improvement of approximately 8.6% in a high-dimensional environment with significant noise interference. Compared to traditional PAM (0.705), it achieves a 24.3% improvement, while K-shape receives a negative score (-0.39) on this dataset due to the mismatch between its time-series shape assumption and the actual pattern. This indicates that LASSO-PAM++ does not represent marginal optimization but rather a substantial lead in high-dimensional clustering tasks. MovingPictures data show that the microbial community in the same location fluctuates significantly on monthly, weekly, and even daily scales. At the same time, most operational taxa do not appear continuously at all time points, indicating that the real population dynamics include both core stable trends and perturbation events.

[0110] Research in high-dimensional clustering within the statistical community emphasizes that sparse clustering strategies, which improve clustering quality through feature selection, have both theoretical and practical advantages when the true cluster differences are driven by a few features. For example, Witten and Tibshirani proposed using an objective function with LASSO penalty within a sparse clustering framework to simultaneously optimize feature selection and cluster assignment, which can more accurately reflect the true cluster structure under high-dimensional conditions. Further theoretical work has also demonstrated that introducing feature selection into high-dimensional mixture models can achieve the minimum clustering error rate with only logarithmic-level dimensionality dependence. LASSO-PAM++ effectively suppresses the interference of noise variables by using L1 regularization for feature selection and then employing PAM's Medoid robust clustering mechanism, thereby achieving a more compact intra-cluster structure and clearer inter-cluster separation, which is the core reason for its high silhouette coefficient.

[0111] Biological Interpretation: Revealing the Dynamic Structure of Microbial Communities Mapping the clustering results of LASSO-PAM++ back to the biological context of the Moving Pictures data allows for the construction of a clear and compelling ecological narrative. The original Moving Pictures data shows that while stable community differences exist between different body parts (such as the gut, mouth, and left and right palms), the microbial communities within a single body part exhibit significant dynamic fluctuations over time. Furthermore, only a few operational taxa appear at high abundance at all time points, suggesting the absence of a truly abundant "core temporal microbiome." Instead, it indicates the presence of persistent but not permanent members and transient species, exhibiting a dual dynamic characteristic of both homeostasis and perturbation.

[0112] LASSO-PAM++ clustering reveals the specific structure of this dual dynamic: approximately 94.5% of the samples are grouped into a dominant cluster, corresponding to the core homeostatic dynamic pattern of the microbial community—the ecological baseline around a common trajectory over time. A few variable clusters (e.g., approximately 4.0% of Cluster 6) reflect a systematic deviation from the core, possibly in response to external disturbances or environmental stimuli. A very small number of clusters (approximately 1.5%) more accurately reflect transient dynamics under acute disturbances or extreme individual behaviors. Within the ecosystem dynamics framework, this discrete structure supports the ecological view of state threshold transitions, suggesting that community states are not continuously and slowly changing, but rather can rapidly transition between different homeostatic states. Compared to clustering methods without feature selection, which may obscure noise from true dynamics, LASSO-PAM++'s feature selection process reduces the interference of noise variables, making the clustering results easier to interpret and more likely to reflect real biological processes. This provides a more robust basis for linking community dynamics with host physiological states or intervention responses.

[0113] This application also provides an electronic device, which includes a processor and a memory. The memory stores computer programs, such as program code implementing the aforementioned sparse microbiome time-series data analysis method. The processor executes the computer program stored in the memory to implement the steps in the aforementioned sparse microbiome time-series data analysis method. The processor can be a central processing unit (CPU), digital signal processor (DSP), application-specific integrated circuit (ASIC), or off-the-shelf programmable gate array (FPGA) or other chip with data processing capabilities. The memory can be a storage medium such as random access memory (RAM), read-only memory (ROM), flash memory, or solid-state drive. The processor and memory can be connected and communicate via a bus. The electronic device may also include necessary components such as communication interfaces and input / output devices.

[0114] This application also provides a computer-readable storage medium storing a computer program thereon. When executed by a processor, the computer program implements the steps of the sparse microbiome time-series data analysis method described in the above embodiments. The computer-readable storage medium can be, but is not limited to, cloud storage servers, USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks, and other media capable of storing program code.

[0115] In the description of this specification, the references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of this application. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0116] The foregoing has shown and described the basic principles, main features, and advantages of this application. Those skilled in the art should understand that this application is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of this application. Various changes and modifications can be made to this application without departing from the spirit and scope thereof, and all such changes and modifications fall within the scope of the claims of this application.

Claims

1. A method for analyzing sparse microbiome time series data, characterized in that, The method includes: Obtaining the time-series data matrix of the microbiome The set number of clusters and regularization parameters The time series data matrix Include Samples of each operational classification unit in The values ​​at each aligned time point Characterizing the time series data matrix The first in Samples of each operational classification unit; In response to the acquired time series data matrix and the number of clusters From the time series data matrix Extract initialization from Each center point set of center points And based on the time series data matrix The time dimension attribute initialization includes Each weight value Non-negative time dimension weight vector ; In response to the acquired time series data matrix The initial set of center points and the non-negative time dimension weight vector The algorithm iteratively solves the problem by minimizing the following objective function, determining whether the parameter changes satisfy a preset convergence condition. If the result indicates that the preset convergence condition is met, the current clustering partition is output. The set of center points and the non-negative time dimension weight vector As the optimal combination of solutions to the objective function; The iterative solution includes: Obtain the samples of the operation classification unit To each of the aforementioned center points The weighted distance is used to classify the samples of each of the operation classification units. Uniquely assigned to the nearest center point The corresponding cluster In the process, the clustering partitions are generated. ; In response to the generated cluster partition From the time series data matrix Extract each of the clusters respectively The local median of the samples of the operation classification unit within the whole at the corresponding time point, and the global median of the samples of the overall operation classification unit at the corresponding time point, are used to generate the inter-cluster sum of squares for each time dimension. In response to the generated inter-cluster sum of squares corresponding to each time dimension and the obtained regularization parameter Generate the updated non-negative time dimension weight vector. ; In response to the updated nonnegative time dimension weight vector and the clustering partition In each of the aforementioned clusters The center point is then determined by re-identifying the actual operational classification unit sample. To update the set of center points .

2. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, The acquisition of time-series data matrix of microbiome ,include: Extract the raw operational taxonomic unit abundance matrix from the acquired microbiome sequencing data. and the corresponding metadata table ; wherein, the original operational classification unit abundance matrix The elements are the absolute counts of each operational taxonomic unit in each sequencing sample, and the metadata table is... Includes a unique sample identifier, gender, and sampling time for each sequencing sample. ; The filtered abundance matrix The absolute count values ​​in the matrix are divided by the total sequencing depth of the corresponding sample to convert them into a relative abundance matrix. ; From the metadata table Extracting the discrete sampling time And by linear scaling the sampling time Mapped to a pre-defined unified time interval, a standardized time vector is generated. ; Based on the relative abundance matrix and standardized time vector Perform linear mixture model spline fitting to generate time series data matrix X.

3. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, The time series data matrix Extract initialization from Each center point set of center points ,include: From the time series data matrix The sample of the first operational classification unit is selected and established as the first centroid and stored in the current set; For each of the unselected operation classification units, the samples Calculate the weighted minimum distance from it to the current set; In response to the calculated weighted minimum distances, a division operation is performed to generate samples for each of the unselected operational classification units. The normalized probability; In response to the generated normalized probability, the sample of the next operational classification unit is extracted, established as the new centroid, and added to the current set. This process is repeated until a set containing the centroid is generated. Each center point The set of center points .

4. The method for analyzing sparse microbiome time series data as described in claim 3, characterized in that, The weighted minimum distance and the normalized probability The calculation formulas are as follows: in, The current set, representing the set generated in the current iteration round, Characterizes the known center point within the current set. The characterization is based on the initialized non-negative time-dimensional weight vector. Calculate the squared distance. Samples representing all unselected operational classification units The sum of the weighted minimum distances.

5. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, The objective function is: in, The characterization includes Cluster Clustering partitioning, Samples characterizing the operational classification unit To its cluster center point The weighted squared distance; In response to the acquired time series data matrix The initial set of center points and the non-negative time dimension weight vector The algorithm iteratively solves the problem by minimizing the following objective function, determining whether the parameter changes satisfy a preset convergence condition. If the result indicates that the preset convergence condition is met, the current clustering partition is output. The set of center points and the non-negative time dimension weight vector This serves as the optimal combination of solutions for solving the objective function.

6. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, The formula for generating the inter-cluster sum of squares corresponding to each time dimension is as follows: in, are positive integers and , Characterizing the first Inter-cluster sum of squares generated by an aligned time dimension Characterizing the cluster Number of samples within, Characterizing the cluster Samples of the operation classification unit within In the The local median at each aligned time point The samples of the operational classification unit that characterize the whole are in the first... The global median at each aligned time point.

7. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, In response to the generated inter-cluster sum of squares corresponding to each time dimension and the obtained regularization parameter Generate the updated non-negative time dimension weight vector. ,include: In response to the generated inter-cluster sum of squares corresponding to each time dimension and the obtained regularization parameter The non-negative time dimension weight vector The update calculation formula is: in, Characterizing the first The updated weight values ​​for each alignment time dimension , and Respectively characterize the first The and the first The inter-cluster sum of squares over each alignment time dimension, Characterize the soft threshold operator.

8. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, In each of the clusters The center point is then determined by re-identifying the actual operational classification unit sample. The formula for the change of center point is: in, Characterizing the cluster The updated center point , and Each represents a group belonging to the same cluster. Samples of two operational classification units within Examples, Characterized by its updated non-negative time dimension weight vector The weighted squared distance value.

9. The method for analyzing sparse microbiome time series data as described in claim 1, characterized in that, Whether the change in the determination parameter satisfies the preset convergence condition includes: Preset tolerance threshold for retrieving storage ; Determine whether the underlying state parameters satisfy the following convergence criterion formula: in, Characterizing the non-negative time dimension weight vector In the current iteration round Compared to the previous iteration round The infinite norm change between them Characterizing the set of center points The change in the number of replacements. Characterizing the clustering partition The change in the number of samples transferred.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the sparse microbiome time series data analysis method as described in any one of claims 1 to 9.

Images