Unsupervised clustering analysis method for molecular subtypes of disease markers

Abstract

The invention provides a disease marker molecular subtype unsupervised clustering analysis method, and belongs to the technical field of disease marker molecular subtypes. Blood samples are collected to obtain tens of thousands of dimensional molecular marker data, and a double-layer game optimization model is constructed for parameter optimization after quality control and preprocessing; the method comprises the following steps: reducing high-dimensional features to a key feature subset by adopting stepped hierarchical sparse coding, solving an optimal grid division scheme by utilizing a Hungary algorithm, constructing a kernel space similarity matrix based on a Gaussian kernel function to redefine distance measurement, and realizing clustering analysis by adopting a spectral clustering algorithm. And performing biological function verification on the clustering result through gene enrichment analysis and pathway analysis, establishing a subtype stability evaluation mechanism to monitor the clustering quality, and finally constructing a subtype specific molecular marker spectrum for clinical diagnosis and prognosis evaluation. The technical problems that distance measurement of a traditional clustering algorithm fails and the clustering effect is sharply reduced due to the fact that the dimension of the molecular feature of the disease marker is extremely high are solved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of disease marker molecule subtypes, and in particular relates to a disease marker molecule subtype unsupervised clustering analysis method. BACKGROUND

[0002] Disease marker molecule subtype analysis is a core technology in the field of precision medicine. Traditional technologies mainly use K-means clustering algorithms based on Euclidean distance, DBSCAN algorithms based on density, and hierarchical clustering algorithms to cluster analyze disease-related molecular markers. By analyzing gene expression profiles, proteomes, and metabolome data, different disease subtype molecular feature patterns are identified. These methods are widely used in clinical application scenarios such as tumor molecular typing, neurodegenerative disease marker screening, cardiovascular disease risk assessment, immune system disease diagnosis, and drug target discovery. Traditional clustering algorithms perform well when dealing with low-dimensional data, but when faced with high-dimensional molecular marker data containing tens of thousands of gene expression levels, protein concentrations, and metabolite contents, they are severely affected by the curse of dimensionality phenomenon. All sample points tend to be equidistantly distributed in high-dimensional space, resulting in a loss of distinguishing ability in similarity calculation based on traditional distance metrics such as Euclidean distance and Manhattan distance. Clustering algorithms cannot effectively identify the inherent structure and patterns in the data. In current disease marker molecule subtype analysis applications, molecular marker data usually contains thousands to tens of thousands of feature dimensions, while the number of samples is relatively small, forming a typical high-dimensional small sample problem. Traditional clustering methods fail in distance measurement under this data structure, and the similarity between all samples tends to be consistent, making it impossible to form meaningful clustering results. The clustering effect decreases sharply with increasing dimensionality, severely restricting the accuracy of molecular subtype identification and the clinical application value. SUMMARY

[0003] Therefore, the application provides a disease marker molecule subtype unsupervised clustering analysis method, which can solve the technical problem of distance measurement failure and sharp decline in clustering effect caused by extremely high molecular feature dimension of disease markers in the prior art.

[0004] The application is implemented in the following manner: the application provides a disease marker molecule subtype unsupervised clustering analysis method, which comprises the following steps: collecting blood samples of disease patients and blood samples of healthy controls, obtaining molecular marker original data comprising gene expression profile data, proteome data and metabolome data through high-throughput sequencing technology, performing quality control detection and nucleic acid extraction and protein separation on the blood samples, establishing a sequencing quality evaluation function to optimize sequencing parameters, performing batch effect correction and preprocessing on the molecular marker original data to obtain a standardized molecular marker data matrix classified into high expression classification, medium expression classification and low expression classification according to expression levels, establishing a double-layer game optimization model to optimize grid division parameters, the upper model takes minimizing clustering error as the target, the lower model takes maximizing feature retention rate as the target, and the two objective functions are associated through a coupling term, performing feature selection based on the standardized molecular marker data matrix by using a ladder hierarchical sparse coding method, constructing a grid optimization problem as a bipartite graph maximum weight matching problem, and solving the optimal grid division parameters by using the Hungarian algorithm, constructing a similarity matrix, and performing unsupervised clustering analysis by using a spectral clustering algorithm, verifying the clustering results in terms of biological functions, establishing a subtype stability evaluation mechanism, and establishing a subtype-specific molecular marker spectrum.

[0005] In the application, the high-throughput sequencing technology refers to a technology for sequencing a large number of DNA or RNA fragments in parallel, which can simultaneously detect tens of thousands to millions of molecular markers in a single experiment.

[0006] In the application, the quality control detection specifically comprises hemoglobin concentration determination, hemolysis degree evaluation and storage time verification, and blood samples meeting the quality standards are selected for subsequent analysis.

[0007] In the application, the nucleic acid extraction and protein separation specifically comprises purifying the extract by using column chromatography and gel electrophoresis to obtain a high-purity molecular marker extract.

[0008] In the application, the column chromatography refers to a chromatographic technology for separating a mixture by using the difference in interaction between a stationary phase and a mobile phase.

[0009] In the application, the gel electrophoresis refers to a technology for separating molecules by using the difference in size and charge under the action of an electric field.

[0010] In the application, the sequencing quality evaluation function specifically comprehensively evaluates the quality of sequencing data, the sequencing depth parameter is adjusted when the evaluation function value is less than a first threshold value, the read length parameter is adjusted when the evaluation function value is less than a second threshold value, and the quality threshold value parameter is improved when the evaluation function value is less than a third threshold value.

[0011] The batch effect correction refers to a data processing method for eliminating systematic differences between different experimental batches, which eliminates systematic errors between different batches through standardization methods to obtain corrected molecular marker data.

[0012] The high expression classification refers to a molecular marker category with an expression level above the upper quartile of the data distribution; the medium expression classification refers to a molecular marker category with an expression level near the median of the data distribution; and the low expression classification refers to a molecular marker category with an expression level below the lower quartile of the data distribution.

[0013] The stepwise hierarchical sparse coding method refers to a coding method for gradually filtering important features from coarse to fine by setting multiple sparsity levels, which gradually filters key features by setting different sparsity thresholds to obtain a feature subset after dimensionality reduction.

[0014] The coupling term refers to a mathematical expression connecting the upper model and the lower model, which realizes the mutual influence between the two models through the product form of the grid size parameter and the retention ratio parameter.

[0015] The Hungarian algorithm is a combinatorial optimization algorithm for solving assignment problems, which finds the maximum weight perfect matching in a bipartite graph by constructing an augmented path. The input includes the grid node weight matrix, the feature similarity matrix, and the constraint condition matrix, and the output is the optimal grid division scheme.

[0016] The similarity matrix is calculated using a Gaussian kernel function to measure the similarity between samples, which is used to quantify the similarity of molecular features between samples.

[0017] The spectral clustering algorithm refers to a clustering method based on graph theory, which realizes clustering by constructing a similarity graph and analyzing its eigenvectors of Laplacian matrix, and determines the number of clusters and cluster centers by eigenvalue decomposition of Laplacian matrix.

[0018] The biological function verification is specifically to evaluate the biological significance of the cluster subtype through gene enrichment analysis and pathway analysis to obtain the classification results of disease marker molecular subtypes.

[0019] The subtype stability evaluation mechanism is specifically to calculate the transition rate of each expression classification, and when the transition rate exceeds a preset threshold, the double-layer game optimization model is returned to optimize the grid division parameters again, and when the transition rate is within the range, the subtype-specific molecular marker spectrum is established.

[0020] The application effectively solves the problem of distance metric failure in high-dimensional molecular marker data by constructing a double-layer game optimization model combined with a ladder hierarchical sparse coding method. The double-layer game model optimizes the feature retention rate in the upper layer and maximizes the feature retention rate in the lower layer, ensuring the data clustering and preserving the key biological feature information in the dimension reduction process. The ladder hierarchical sparse coding method sets multiple sparsity thresholds to filter step by step, compressing the tens of thousands of original feature space to hundreds of effective feature space, avoiding information loss caused by one-time dimension reduction. The application uses the Hungarian algorithm to solve the optimal grid partition parameters, ensuring the global optimality of the feature selection process, and then uses the Gaussian kernel function to map the reduced features to the kernel space, redefining the similarity measurement between samples in the kernel space and overcoming the problem of distance metric failure in the original high-dimensional space. The spectral clustering algorithm converts the clustering problem into an optimal partition problem in graph theory by analyzing the Laplacian matrix eigenvectors of the similarity matrix, avoiding the direct dependence of traditional methods on distance measurement. The subtype stability evaluation mechanism established by the application monitors the expression classification transition rate to evaluate the stability of the dimension reduction and clustering process in real time, automatically triggers parameter re-optimization when the clustering effect is detected to be declining, and forms an adaptive quality control system, solving the technical problems of distance metric failure and sharp decline in clustering effect caused by the high dimension of disease marker molecular features in the prior art. BRIEF DESCRIPTION OF DRAWINGS

[0021] Figure 1 A flowchart of the method of the application.

[0022] Figure 2 A convergence process diagram of the double-layer game optimization model in Example 2.

[0023] Figure 3 A comparison chart of the number of molecular subtype-specific markers in Example 2.

[0024] Figure 4 A sequencing quality evaluation function change chart in Example 2.

[0025] Figure 5 A sample similarity distribution scatter plot in Example 2.

[0026] Figure 6 A feature screening process chart of the ladder hierarchical sparse coding in Example 2.

[0027] Figure 7 A clustering result stability evaluation chart in Example 2.

[0028] Figure 8 A feature value distribution and cluster number determination chart in Example 2. DETAILED DESCRIPTION

[0029] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application.

[0030] As shown in Figure 1 is a flow chart of a method of unsupervised clustering analysis of disease marker molecular subtypes provided by the present application, the method comprises the following steps:

[0031] S01, collecting blood samples of disease patients and blood samples of a healthy control group, centrifugally separating the blood samples to obtain serum and plasma components, and obtaining molecular marker original data through high-throughput sequencing technology, wherein the molecular marker original data comprises gene expression profile data, proteome data and metabolome data;

[0032] S02, performing quality control detection on the blood samples, including hemoglobin concentration determination, hemolysis degree evaluation and storage time verification, and screening blood samples meeting the quality standards for subsequent analysis;

[0033] S03, performing nucleic acid extraction and protein separation on the blood samples meeting the quality standards, purifying the extracts by column chromatography and gel electrophoresis to obtain high-purity molecular marker extracts;

[0034] S04, establishing a sequencing quality evaluation function to optimize the sequencing parameters of the high-purity molecular marker extracts, adjusting the sequencing depth parameter when the evaluation function value is less than a first threshold value, adjusting the read length parameter when the evaluation function value is less than a second threshold value, and improving the quality threshold parameter when the evaluation function value is less than a third threshold value;

[0035] S05, performing batch effect correction and technical repetition processing on the molecular marker original data, eliminating systematic errors between different batches by a standardization method to obtain corrected molecular marker data;

[0036] S06, preprocessing the corrected molecular marker data, including data standardization, missing value filling and outlier detection, to obtain a standardized molecular marker data matrix, wherein the standardized molecular marker data matrix is divided into high expression classification, medium expression classification and low expression classification according to expression level;

[0037] S07, establishing a double-layer game optimization model to optimize grid division parameters, the upper model taking minimizing clustering error as the target, the lower model taking maximizing feature retention rate as the target, and the two objective functions being associated through a coupling term;

[0038] S08, based on the standardized molecular marker data matrix, using a ladder hierarchical sparse coding method to select features, and screening key features step by step by setting different sparsity threshold values to obtain a feature subset after dimension reduction.

[0039] S09, constructing the grid optimization problem as a maximum weight bipartite matching problem, solving the optimal grid partition parameters by the Hungarian algorithm, and obtaining the optimized grid partition scheme;

[0040] S10, constructing a similarity matrix for the dimensionality-reduced feature subset, and calculating the similarity measure between samples by using a Gaussian kernel function, the similarity measure being used to quantify the similarity degree of molecular features between samples;

[0041] S11, performing unsupervised clustering analysis based on the similarity matrix by using a spectral clustering algorithm, and determining the number of clusters and cluster centers by eigenvalue decomposition of a Laplacian matrix;

[0042] S12, verifying the clustering results biologically, evaluating the biological significance of the clustered subtypes by gene enrichment analysis and pathway analysis, and obtaining the classification results of disease marker molecular subtypes;

[0043] S13, establishing a subtype stability evaluation mechanism, calculating the transition rate of each expression classification, returning to step S07 to re-optimize the grid partition parameters when the transition rate exceeds a preset threshold, and entering step S14 when the transition rate is within an acceptable range;

[0044] S14, establishing a subtype-specific molecular marker spectrum, sorting and screening the feature molecules of each subtype, and constructing a molecular marker combination for clinical diagnosis and prognosis evaluation.

[0045] The high-throughput sequencing technology refers to a technology for simultaneously detecting tens of thousands to millions of molecular markers in a single experiment by parallel sequencing of a large number of DNA or RNA fragments.

[0046] The column chromatography refers to a chromatographic technology for separating a mixture by using the difference in interaction between a stationary phase and a mobile phase.

[0047] The gel electrophoresis method refers to a technical method for separating by using the difference in molecular size and charge under the action of an electric field.

[0048] The mathematical expression of the sequencing quality evaluation function is F seq =α·D depth +β·L read +γ·Q score +δ·C coverage , which is used for comprehensive evaluation of sequencing data quality, wherein α, β, γ, and δ are weight coefficients.

[0049] The sequencing depth parameter refers to the average number of times each target region is sequenced.

[0050] The read length parameter refers to the length of the nucleotide sequence read each time.

[0051] Quality score parameter refers to the quantitative indicator of sequencing accuracy.

[0052] Coverage parameter refers to the percentage of target genomic regions covered by sequencing.

[0053] Batch effect correction refers to a data processing method to eliminate systematic differences between different experimental batches.

[0054] High expression classification refers to the class of molecular markers with expression levels above the upper quartile of the data distribution.

[0055] Medium expression classification refers to the class of molecular markers with expression levels around the median of the data distribution.

[0056] Low expression classification refers to the class of molecular markers with expression levels below the lower quartile of the data distribution.

[0057] The mathematical expression of the objective function of the upper model in the double-layer game optimization model is where x i is the coordinate of the cluster center, c j is the weight of the grid node, θ i is the angle parameter, d i is the distance parameter, and λ, μ, v, γ are adjustment coefficients.

[0058] The mathematical expression of the objective function of the lower model in the double-layer game optimization model is where u k is the retention ratio parameter, v k is the redundancy parameter, φ k is the phase parameter, r k is the radius parameter, and ρ, σ, τ, ω are adjustment coefficients.

[0059] The coupling term refers to the mathematical expression connecting the upper model and the lower model, which realizes the mutual influence between the two models through the product form of the grid size parameter and the retention ratio parameter.

[0060] Contribution rate refers to the quantitative indicator of the influence degree of each expression classification on the final clustering result.

[0061] Transition rate refers to the frequency of molecular markers migrating between different expression classifications.

[0062] Ladder hierarchical sparse coding method refers to a coding method that sets multiple sparsity levels to gradually filter important features from coarse to fine, effectively solving the feature redundancy problem in high-dimensional data.

[0063] Hungarian algorithm is a combinatorial optimization algorithm for solving assignment problems, which finds the maximum weight perfect matching in a bipartite graph by constructing an augmenting path. The input includes the grid node weight matrix, the feature similarity matrix and the constraint condition matrix. The output is the optimal grid partition scheme.

[0064] Gaussian kernel function is a kind of radial basis function used to calculate the nonlinear similarity between samples, which is the product of exponential function and squared Euclidean distance.

[0065] Spectral clustering algorithm is a clustering method based on graph theory, which realizes clustering by constructing a similarity graph and analyzing the eigenvectors of its Laplacian matrix.

[0066] Gene enrichment analysis is an analysis method that evaluates whether a gene set is significantly enriched in a biological process through statistical methods.

[0067] Pathway analysis is an analysis method that studies the interaction relationship of genes or proteins in biological pathways.

[0068] Subtype-specific molecular marker profile refers to the unique combination of molecular characteristics of each disease subtype, which is used to distinguish the set of molecular markers of different subtypes.

[0069] The specific embodiments of the above steps are described in detail as follows.

[0070] The specific embodiment of step S01 is to realize the standardized collection of raw data through systematic blood sample collection and multi-omics data acquisition. First, the venous blood samples of disease patients and healthy controls are collected by using vacuum blood collection tubes, and the blood collection amount is controlled at 10-15 milliliters. The blood collection time is selected at 8-10 am on an empty stomach to reduce the influence of physiological rhythm on the expression of molecular markers. Then the blood samples are left at room temperature for 30 minutes to coagulate naturally, and then centrifuged at 3000 revolutions per minute at 4℃ for 15 minutes to separate serum and plasma components. Serum is used for proteomics analysis, and plasma is used for metabolomics detection. The extracted nucleic acids are subjected to whole genome sequencing by next-generation sequencing technology to obtain gene expression profile data, and liquid chromatography mass spectrometry is used to detect proteomic data, and gas chromatography mass spectrometry is used to obtain metabolomic data. The purpose of this step is to establish a comprehensive molecular marker data set containing three levels of genes, proteins and metabolites, and to provide complete biological information basis for subsequent analysis.

[0071] The specific implementation of step S02 is to establish a multi-dimensional quality control system to ensure the consistency and reliability of sample quality. Hemoglobin concentration determination uses spectrophotometry. After diluting the blood sample 100 times, the absorbance is measured at a wavelength of 540 nanometers. The normal range should be controlled within 120-160 grams per liter. Samples exceeding this range will be excluded. Hemolysis assessment is achieved by detecting the free hemoglobin content in plasma. Benzidine method is used for detection. When the free hemoglobin concentration exceeds 0.2 grams per liter, it is determined as a hemolytic sample and is rejected. Storage time verification requires that the time from blood collection to processing should not exceed 2 hours. The storage temperature of frozen samples should be kept at minus 80 degrees Celsius, and the storage time should not exceed 6 months. Quality control also includes detecting the integrity of the sample DNA. The degree of DNA degradation is evaluated by agarose gel electrophoresis. The intact DNA should present a clear high molecular weight band. This step excludes low-quality samples that may affect the accuracy of subsequent analysis results through strict quality control standards.

[0072] The specific implementation of step S03 is to obtain high-purity molecular marker extracts using multiple separation and purification techniques. Nucleic acid extraction uses phenol chloroform extraction method. First, proteinase K is used to digest proteins at 56°C for 2 hours, then an equal volume of phenol chloroform isopropyl alcohol mixture is added to remove proteins. After multiple extractions, nucleic acids are precipitated with isopropyl alcohol. Finally, high-quality DNA and RNA with a purity ratio of 1.8-2.0 are obtained. Protein separation uses ammonium sulfate fractionation precipitation method combined with ion exchange chromatography technology. First, different types of proteins are precipitated step by step by different concentrations of ammonium sulfate solution, then further separated and purified by anion exchange column. Column chromatography uses gel filtration chromatography principle to separate according to the difference in molecular size. The appropriate pore size of the gel medium is selected to ensure the effective separation of the target molecules. Gel electrophoresis uses the molecular sieve effect of polyacrylamide gel to separate according to the molecular weight under the action of electric field. The separation precision is controlled by adjusting the gel concentration. The purpose of this step is to obtain molecular marker extracts with a purity higher than 95%, providing a reliable sample basis for subsequent high-throughput detection.

[0073] The specific implementation of step S04 is to construct a comprehensive evaluation function to realize dynamic optimization adjustment of sequencing parameters. The sequencing quality evaluation function comprehensively considers four key parameters of sequencing depth, read length, quality score and coverage, and calculates the comprehensive evaluation value by weighted summation. The first threshold is set to 0.75, and when the evaluation function value is lower than this threshold, it indicates that the sequencing depth is insufficient, and the sequencing depth needs to be increased from the initial 30 times to 50 times or higher. The second threshold is set to 0.65, and when the evaluation function value is lower than this threshold, it indicates that the read length parameter needs to be adjusted, and the 150 base pair read length is extended to 250 base pairs to obtain more complete sequence information. The third threshold is set to 0.55, and when the evaluation function value is lower than this threshold, the quality threshold parameter needs to be improved, and the Phred quality score is increased from 20 to 30 to ensure that the sequence accuracy is above 99.9%. The weight coefficient is adjusted according to different molecular types, and the sequencing depth weight is higher in gene expression profile data, and the coverage weight is larger in proteome data. This step ensures that high-quality multi-omics data is obtained through dynamic parameter optimization, and provides a reliable data foundation for subsequent analysis.

[0074] The specific implementation of step S05 is to eliminate systematic errors and technical variations by using a multi-level correction strategy. Batch effect correction uses the ComBat algorithm, which adjusts the location and scale parameters between different batches based on the empirical Bayes method, and identifies and removes batch-related systematic bias by establishing a linear mixed effect model. Technical repeatability processing uses analysis of variance method to evaluate the consistency between repeated samples, and calculates the intraclass correlation coefficient. When the correlation coefficient is less than 0.9, the experiment is re-performed. The standardization method uses quantile standardization technology to adjust the expression distribution of all samples to the same quantile distribution, eliminating technical differences between samples. For gene expression data, the DESeq2 standardization method is used to correct the library size difference by calculating the scale factor. The proteome data uses total protein standardization method, which normalizes the total protein concentration as the reference. The metabolome data uses internal standard normalization method, which adds internal standard compounds with known concentration to correct the variation in the detection process. This step ensures that the data generated by different batches and technical platforms is comparable, laying a foundation for subsequent integrated analysis.

[0075] The specific implementation of step S06 is to establish a multi-dimensional data preprocessing process to further improve the data quality. The data standardization adopts Z-score standardization method to convert the value of each feature to a standard normal distribution with mean of 0 and standard deviation of 1, eliminating the influence of dimension difference between different features. The missing value filling adopts K-nearest neighbor interpolation algorithm, which selects the nearest K neighbor samples according to the similarity between samples, estimates the missing value by weighted average method, and the K value is set to 5, and the weight is determined according to the inverse of Euclidean distance. The outlier detection adopts the method based on quartile range, calculates the first quartile and third quartile of each feature, and judges the values exceeding 1.5 times of the quartile range as outliers and processes them. The expression level classification determines the boundary point according to the statistical characteristics of data distribution, and the high expression classification corresponds to the values above the 75th percentile, the medium expression classification corresponds to the values between the 25th percentile and the 75th percentile, and the low expression classification corresponds to the values below the 25th percentile. This classification method can effectively capture the dynamic range of molecular marker expression and provide meaningful feature division for subsequent clustering analysis.

[0076] The specific implementation of step S07 is to construct a double-layer game optimization framework to realize intelligent optimization of grid division parameters. The upper model takes minimizing clustering error as the goal, and optimizes the clustering center position, grid node weight, angle parameter and distance parameter to improve the clustering quality. The objective function considers the combination of Euclidean distance square, logarithmic function, sine function and exponential function, which can effectively handle nonlinear and multi-modal optimization problems. The lower model takes maximizing feature retention rate as the goal, and adjusts the retention ratio parameter, redundancy parameter, phase parameter and radius parameter to effectively retain important features. The objective function adopts the combination of linear term, quadratic term, cosine function and logarithmic function. The coupling term establishes the interaction relationship between the upper and lower models through the product form of grid size parameter and retention ratio parameter, ensuring the coordinated development of the two optimization objectives. The game process adopts the alternating optimization strategy, which first optimizes the upper model by fixing the lower variable, then optimizes the lower model by fixing the upper variable, and iterates until convergence. This double-layer optimization framework can consider the dual requirements of clustering accuracy and feature retention, and realize more balanced and effective parameter optimization.

[0077] The specific implementation of step S08 is to realize hierarchical screening of high-dimensional features by using a stepwise hierarchical sparse coding technique. The method sets the sparsity threshold to multiple increasing levels, starting from 0.9 and decreasing by 0.1 to 0.1, forming 9 sparsity levels. Each sparsity level uses the L1 regularization sparse coding algorithm to obtain sparse representation coefficients by solving a constrained optimization problem. At the coarsest sparsity level 0.9, preliminary feature screening is performed, retaining the top 10% features with the largest absolute value of coefficients, and then further screening the retained features at the next sparsity level 0.8. During the stepwise screening process, the stability of the feature subset is evaluated using the cross-validation method, and features with a reproduction rate higher than 80% in multiple resampling are selected as candidate features. The feature importance score is evaluated comprehensively based on the sparsity coefficient size, frequency of occurrence and biological significance. This stepwise screening strategy can effectively avoid the loss of important information caused by one-time feature selection, while controlling the feature dimension within a manageable range, providing the most representative feature subset for subsequent clustering analysis.

[0078] The specific implementation of step S09 is to convert the grid optimization problem into a bipartite graph maximum weight matching problem and solve it using the Hungarian algorithm. First, a bipartite graph structure is constructed, with grid nodes as one side of vertices and feature subsets as the other side of vertices. The weight of the edge is calculated according to the grid node weight matrix, feature similarity matrix and constraint condition matrix. The grid node weight matrix reflects the contribution of different grid positions to the clustering effect, and the weight distribution is determined by statistical analysis of historical clustering results. The feature similarity matrix calculates the Pearson correlation coefficient between features, and features with high similarity are assigned to adjacent grid nodes to maintain local consistency. The constraint condition matrix contains the boundary conditions and balance requirements of grid division, ensuring that each grid region contains an appropriate number of features. The Hungarian algorithm gradually expands the matching scale by constructing an augmented path until the maximum weight perfect matching is found. The time complexity of the algorithm is O(n 3 ), where n is the number of nodes. This method can guarantee the global optimal grid division scheme under various constraint conditions, providing an optimized data organization structure for subsequent clustering analysis.

[0079] The specific implementation of step S10 is to construct a similarity matrix based on a Gaussian kernel function to quantify the similarity of molecular features between samples. The Gaussian kernel function adopts a radial basis function form, and the core idea is to take the square of the Euclidean distance as the negative exponent of the exponential function, and to control the similarity decay rate by adjusting the bandwidth parameter. The bandwidth parameter is determined by empirical rules, and is set to the inverse of the square root of the feature dimension. For a 1000-dimensional feature space, the bandwidth parameter is about 0.032. In the similarity calculation process, the Euclidean distance between sample pairs is first calculated, and then the similarity measurement value is obtained by substituting the Gaussian kernel function formula. To improve the calculation efficiency, the fast kernel method is used, and the Nystrom approximation technique is used to reduce the computational complexity. The similarity matrix has symmetry and positive definiteness, and the diagonal elements are 1, indicating that the sample is completely similar to itself, and the non-diagonal elements are between 0 and 1, indicating the similarity between different samples. This similarity measurement can effectively capture the nonlinear similarity relationship between samples in high-dimensional space, and provide high-quality input data for the spectral clustering algorithm.

[0080] The specific implementation of step S11 is to use a spectral clustering algorithm to realize unsupervised clustering analysis based on the feature decomposition of the similarity matrix. First, a similarity graph is constructed according to the similarity matrix, and a K-nearest neighbor graph construction method is used to select the K neighbors with the highest similarity for each sample to establish a connection, and the value of K is set to the square root of the number of samples. Then, the Laplacian matrix of the graph is calculated, and the standardized Laplacian matrix form is used to avoid the singularity problem of the degree matrix. The Laplacian matrix is decomposed, and the eigenvectors corresponding to the smallest number of eigenvalues are calculated, and the number of eigenvectors is equal to the expected number of clusters. The number of clusters is determined by the feature gap heuristic method, and the position where a significant jump occurs in the eigenvalue sequence is selected as the optimal cluster number. The matrix composed of the eigenvectors is clustered by K-means clustering by row to obtain the final clustering result. The spectral clustering algorithm can handle arbitrary shaped cluster structures, and is particularly suitable for handling complex clustering patterns in high-dimensional molecular data. This method realizes effective identification of molecular subtypes by combining graph theory and linear algebra.

[0081] The specific implementation of step S12 is to evaluate the biological significance and clinical relevance of the clustered subtypes through multi-level biological function verification. Gene enrichment analysis uses the hypergeometric test method to compare the differentially expressed genes in each clustered subtype with known biological function gene sets, calculate the significance level of enrichment, and gene sets with P value less than 0.05 are considered to be significantly enriched. The functional gene set includes three levels of biological process, molecular function and cellular component in the gene ontology database, and metabolic pathway data in the Kyoto Encyclopedia of Genes and Genomes. Pathway analysis uses a topological method to consider the position and interaction relationship of genes in the pathway, and calculates the pathway activity score to evaluate the difference in pathway activity between different subtypes. Protein-protein interaction network analysis identifies key regulatory molecules and signal pathway hubs by constructing a protein-protein interaction graph. Metabolic pathway analysis focuses on the upstream and downstream relationships of metabolites and enzyme-catalyzed reactions, and identifies metabolic characteristic differences between different subtypes. This comprehensive verification strategy ensures that the identified molecular subtypes have clear biological basis and potential clinical value.

[0082] The specific implementation of step S13 is to establish a dynamic stability evaluation mechanism to monitor the reliability and consistency of the subtype classification. The transition rate is calculated by comparing the classification results at consecutive time points or on different sample sets, and the frequency of migration of molecular markers between high expression, medium expression and low expression classifications is calculated. The calculation formula of the transition rate is the number of classification changes divided by the total number of molecules, and the preset threshold is set to 15%. This threshold is determined based on statistical analysis of multiple independent experiments, and when the transition rate is less than 15%, the classification result is considered stable and reliable. The contribution rate evaluation calculates the influence degree of each expression classification on the clustering result, uses leave-one-out cross-validation, removes each expression classification in turn and re-clusters, and compares the degree of change in the clustering result. Stability evaluation also includes resampling verification, which generates multiple resampling data sets through the bootstrap method to evaluate the reproducibility of the clustering result. When the transition rate exceeds the preset threshold, it means that the current grid division parameters are not optimized enough, and it is necessary to return to step S07 to re-optimize the parameters, adjust the weight coefficients and constraint conditions in the game model. This dynamic monitoring mechanism ensures the stability and repeatability of the analysis results.

[0083] The embodiment of step S14 is to construct subtype-specific molecular marker profiles for clinical diagnosis and prognosis evaluation applications. The importance of the characteristic molecules of each subtype is ranked, taking into account three dimensions of expression fold difference, statistical significance and biological function correlation. The expression fold difference uses the log-transformed fold change value, and the molecular marker with significant difference is required to have a fold change greater than 2 times. The statistical significance uses the P value corrected by multiple hypothesis testing, and the Benjamini-Hochberg method is used to control the false discovery rate below 5%. The biological function correlation is evaluated by literature mining and database annotation to assess the strength of the known association between the molecular marker and the disease. The construction of the molecular marker combination uses the stepwise regression method, starting from the most important molecule and adding one by one until the diagnostic performance is no longer significantly improved. The diagnostic performance is evaluated by the area under the receiver operating characteristic curve, and the area under the curve of the diagnostic marker combination is required to be greater than 0.85. The selection of prognostic evaluation markers focuses on molecules related to survival time and treatment response, and the Cox proportional hazards model is used to evaluate the prognostic value. The molecular marker profile provides important molecular basis for precise diagnosis and individualized treatment of different disease subtypes.

[0084] It should be noted that the key technical ideas of the present application mainly reflect in four aspects of double-layer game optimization model, ladder hierarchical sparse coding method, multi-omics data integration analysis and dynamic stability evaluation mechanism.

[0085] The double-layer game optimization model optimizes the clustering accuracy and feature retention rate, two mutually competitive objectives, by constructing an optimization framework for the interaction between the upper and lower layers. Compared with traditional single-objective optimization methods, this model can maximize the retention of important biological features while ensuring clustering quality, avoiding the loss of important information in the feature selection process. The upper model focuses on the accuracy of the clustering results, and the lower model focuses on the integrity of the feature information. Through the coupling term, a coordination mechanism is established between the two, achieving an effective balance of multi-objective optimization. This game theory method shows stronger robustness and adaptability in dealing with complex optimization problems of high-dimensional biological data.

[0086] The ladder hierarchical sparse coding method realizes step-by-step feature selection from coarse to fine by setting multiple decreasing sparsity thresholds. Traditional one-time feature selection methods are easily affected by noise interference and local optimal traps, while the ladder method can better capture the hierarchical relationships and dependencies between features through a multi-level filtering mechanism. This method performs independent optimization at each sparsity level and then passes the results to the next level, forming a feature selection strategy from coarse to fine, effectively improving the stability and biological interpretability of feature selection.

[0087] The multi-omics data integration analysis systematically integrates gene expression, proteome and metabolome data to construct a more comprehensive and accurate molecular marker atlas. Traditional methods usually only focus on single omics level information, which is difficult to fully reflect the molecular mechanism of the disease. The application realizes the effective fusion of different omics data through a unified data processing and analysis framework, can simultaneously identify disease-related molecular changes from multiple biological levels, and improves the accuracy and biological significance of subtype identification.

[0088] The dynamic stability evaluation mechanism establishes a self-adaptive quality control system by real-time monitoring of the stability of the classification results. This mechanism can timely find the unstable factors that may occur in the analysis process, trigger the parameter re-optimization process through the quantitative evaluation of the conversion rate and contribution rate, and ensure the reliability of the analysis results. Compared with the static quality control method, the dynamic evaluation mechanism has stronger self-adaptive ability and error correction ability.

[0089] The synergistic effect of these four key technical ideas forms a complete intelligent molecular subtype identification system. The double-layer game model provides an optimized parameter configuration for the entire analysis process, the ladder grading coding ensures the quality of feature selection, the multi-omics integration provides comprehensive molecular information, and the dynamic evaluation guarantees the stability of the results. This multi-technology collaborative approach significantly improves the accuracy, stability and biological interpretability of molecular subtype identification, providing a more reliable molecular basis for precision medicine and individualized treatment. Compared with existing technologies, the application can handle more complex high-dimensional data, identify more accurate molecular subtypes, and provide more stable analysis results, which has important clinical application value in disease diagnosis and treatment selection.

[0090] Specifically, the principle of the present application is that the present application can solve the technical problem that the dimension of disease marker molecules is extremely high, resulting in invalid distance measurement of traditional clustering algorithms and sharp decline of clustering effect, and the technical principle is to overcome the curse of dimensionality in high-dimensional space through the synergistic effect of multi-level dimension reduction optimization and kernel space mapping. The design principle of the double-layer game optimization model is based on multi-objective optimization theory. The upper model ensures the clusterability of the data after dimension reduction by minimizing the clustering error function, and the lower model ensures that important biological information is not lost by maximizing the feature retention rate function. The two objective functions are mathematically related through a coupling term, forming a set of Pareto optimal solutions. This design principle ensures that the key molecular feature information is not lost while solving the high-dimensional problem. The core principle of the ladder hierarchical sparse coding method is to use the sparse induction characteristics of L1 regularization. By setting a decreasing sparse threshold sequence, features are gradually filtered from coarse to fine. First, remove redundant features irrelevant to clustering goals, and then further refine and filter the remaining features. This hierarchical strategy effectively avoids the problem of over-dimension reduction or insufficient dimension reduction that may occur in traditional dimension reduction methods. Meanwhile, the Hungarian algorithm converts the feature selection problem into a combinatorial optimization problem by constructing a bipartite graph matching model of feature importance, ensuring global optimal selection of feature subsets. The mathematical principle of Gaussian kernel function mapping is based on the theory of reproducing kernel Hilbert space. By nonlinear transformation, the original high-dimensional feature space is mapped to an infinite-dimensional kernel space. In the kernel space, originally linearly inseparable high-dimensional data becomes linearly separable. The radial basis characteristics of the kernel function make the inner product value of similar samples in the kernel space larger, and the inner product value of dissimilar samples smaller, re-establishing an effective similarity measurement standard. The spectral clustering algorithm is based on graph theory and matrix analysis theory. By constructing a sample similarity graph and calculating its Laplacian matrix eigenvectors, the clustering problem is converted into an optimal graph cut problem, avoiding direct distance calculation in high-dimensional space. Instead, clustering is achieved by analyzing the connectivity and segmentation cost of the graph. This indirect clustering strategy effectively avoids the fundamental problem of distance measurement failure in high-dimensional space. The various components of the technical solution are organically combined through the complementarity of theoretical foundations and algorithmic processes, forming a complete high-dimensional data clustering solution.

[0091] A specific embodiment 1 of the present application is provided below, and the specific implementation of each step in embodiment 1 is described in detail as follows.

[0092] In this embodiment, the specific implementation of steps S01-S03 is the same as described above, and will not be described in detail here.

[0093] The specific implementation of step S04 is to establish a comprehensive evaluation function to dynamically optimize and adjust the sequencing parameters. The mathematical expression of the sequencing quality evaluation function is:

[0094] F seq = α · Ddepth + β · L read + γ · Q score + δ · C coverage ;

[0095] In the formula, F seq is a sequencing quality comprehensive evaluation value, with a value range of 0-1; D depth is a sequencing depth parameter, representing the average number of times each target region is sequenced, obtained by calculating the number of reads aligned to a certain region of the reference genome; L read is a read length parameter, representing the length of the nucleotide sequence read each time, with a unit of base pairs; Q score is a quality score parameter, representing a quantitative indicator of sequencing accuracy, obtained using the Phred quality score calculation method; C coverage is a coverage parameter, representing the percentage of the target genomic region covered by sequencing; α, β, γ, δ are weight coefficients, respectively taking values of 0.3, 0.25, 0.25, and 0.2. Among them, the parameter acquisition method is: D depth is calculated by counting the number of read coverage times for each genomic site in the alignment file, and the calculation formula is total read number divided by target region length; L read is directly set by the sequencing instrument parameter, and the commonly used value is 150 or 250 base pairs; Q score is calculated through the quality file output by the sequencing instrument, using the formula -10log 10 (P error ), wherein P error is the base calling error probability; C coverage is obtained by calculating the proportion of the number of genomic sites covered by at least one read to the total number of target sites. When F seq is less than the first threshold value 0.75, adjust D depth from 30 times to 50 times; when F seq is less than the second threshold value 0.65, adjust L read from 150 base pairs to 250 base pairs; when F seq is less than the third threshold value 0.55, increase Q score threshold from 20 to 30.

[0096] The specific implementation of step S05 is the same as the foregoing, and will not be described in detail here.

[0097] The specific implementation of step S06 is to establish a standardized molecular marker data matrix and perform expression level classification. The data standardization adopts the Z-score standardization method, and the calculation formula is:

[0098]

[0099] In the formula, Zij is the standardized data matrix element; X ij is the expression value of the jth feature of the ith sample in the original data matrix; μ j is the mean of the jth feature in all samples; σ j is the standard deviation of the jth feature in all samples. The expression level classification adopts the quantile division method, the high expression classification threshold T high = Q 75 , the medium expression classification threshold range is Q 25 ≤ T mid ≤ Q 75 , and the low expression classification threshold T low = Q 25 , where Q 25 and Q 75 are the 25th and 75th percentiles, respectively.

[0100] The specific implementation of step S07 is to construct a double-layer game optimization model to optimize the grid division parameters. The objective function of the upper-layer model is:

[0101]

[0102] In the formula, x i is the coordinate vector of the ith clustering center, with a dimension of the feature space dimension; c i is the weight of the ith grid node, obtained by statistical analysis of historical clustering effects; y i is the number of samples in the ith cluster; θ i is an angle parameter, representing the angle between the clustering center and the coordinate axis, with a value range of 0-2π; d i is a distance parameter, representing the Euclidean distance from the clustering center to the data center; λ, μv, and η are adjustment coefficients, with values of 0.1, 0.05, 0.08, and 1.5, respectively; and n is the number of clusters. The objective function of the lower-layer model is:

[0103]

[0104] In the formula, u k is the retention ratio parameter of the kth feature, with a value range of 0-1; v k is the redundancy parameter of the kth feature, obtained by calculating the Pearson correlation coefficient between features; φ k is a phase parameter, representing the phase information of the feature in the frequency domain space; r k is a radius parameter, representing the distribution radius of the feature in the high-dimensional space; ρξ, τ, and ω are adjustment coefficients, with values of 0.6, 0.2, 0.1, and 0.1, respectively. The coupling term expression is:

[0105]

[0106] In the formula, H couple The coupling strength is w. ik is the grid size parameter, representing the size of the i-th grid cell; m is the total number of features.

[0107] The specific implementation method of step S08 is the same as described above, and will not be repeated in detail here.

[0108] The specific implementation method of step S09 is the same as described above, and will not be repeated in detail here.

[0109] The specific implementation of step S10 involves constructing a similarity matrix based on a Gaussian kernel function. The mathematical expression for the Gaussian kernel function is:

[0110]

[0111] In the formula, K(x) p x q ) is the sample x p and sample x q The similarity measure between them ranges from 0 to 1; x p and x q These are the feature vectors of the p-th and q-th samples, respectively; ||x p -x q || 2 σ is the squared Euclidean distance between two samples; kernel The kernel function bandwidth parameter takes the value of Where D feature The feature dimension is denoted as . The elements of the similarity matrix S are calculated as follows:

[0112] S pq =K(x) p x q );

[0113] In the formula, S pq Let be the element in the p-th row and q-th column of the similarity matrix; the similarity matrix S is an N×N symmetric positive definite matrix, where N is the total number of samples.

[0114] The specific implementation of step S11 involves using a spectral clustering algorithm based on the eigenvalue decomposition of the Laplacian matrix for cluster analysis. The formula for calculating the degree matrix M is:

[0115]

[0116] In the formula, M ii Let L be the diagonal element of the degree matrix, representing the sum of similarities between the i-th sample and all other samples; M is an N×N diagonal matrix. The normalized Laplace matrix L... norm The calculation formula is:

[0117] L norm = M -1 / 2 (M-S)M -1 / 2 ;

[0118] In the formula, L norm is a normalized Laplacian matrix; M -1 / 2 is a negative two-division power of the degree matrix M. The eigenvalue decomposition process is to solve a generalized eigenvalue problem:

[0119] L norm e j = β j e j ;

[0120] In the formula, e j is the jth eigenvector; β j is the corresponding eigenvalue, arranged in ascending order; the number of clusters h is determined by the eigenvalue gap heuristic method, and the value of h is selected to make β h+1 - β h maximal.

[0121] The specific implementation of step S12 is the same as the foregoing, and will not be described in detail here.

[0122] The specific implementation of step S13 is to establish a subtype stability evaluation mechanism. The conversion rate calculation formula is:

[0123]

[0124] In the formula, R trans is the conversion rate; N changed is the number of molecular markers that have changed in expression classification; N total is the total number of molecular markers. The contribution rate calculation formula is:

[0125]

[0126] In the formula, C contrib is the contribution rate; ΔS cluster is the degree of change in clustering results after removing a certain expression classification; ΔS total is the total clustering change degree. When R trans > 0.15, return to step S07 to re-optimize.

[0127] The specific implementation of step S14 is the same as the foregoing, and will not be described in detail here.

[0128] It should be noted that in the present embodiment, the sequencing quality evaluation function F seqThe function adopts the form of multi-parameter weighted summation, and realizes the comprehensive evaluation of sequencing data quality by comprehensively considering four key indicators: sequencing depth, read length, quality score and coverage. The core principle of the function is based on a linear weighted model, and each parameter represents an important quality dimension in the sequencing process, and the importance of each parameter is adjusted by the weight coefficient. Compared with the traditional single indicator evaluation method, the comprehensive evaluation function can more accurately reflect the overall quality of the sequencing data, avoid the quality judgment deviation caused by the abnormality of a single indicator, and thus provide a more reliable data basis for the subsequent molecular marker analysis.

[0129] The upper objective function of the double-layer game optimization model combines four components: distance term, logarithmic term, trigonometric function term and exponential term. Each component undertakes different optimization tasks. The distance term ||x i -c j || 2 is responsible for minimizing the dispersion within the cluster, the logarithmic term log(y i is used to balance the sample distribution of different clusters, the trigonometric function term sin(θ i introduces periodic constraints to avoid local optimum, and the exponential term exp(-γ·d i provides distance decay effect. The lower objective function adopts the maximization strategy, the linear term u k directly optimizes the feature retention ratio, the quadratic term controls the nonlinear growth of redundancy, the cosine term cos(φ k handles the phase relationship between features, and the logarithmic term ln(r k realizes the logarithmic growth of radius parameter. The double-layer structure realizes the dynamic balance between clustering accuracy and feature retention rate through the interaction of upper and lower objective functions, compared with the traditional single objective optimization method, it can maximize the retention of important biological features while ensuring the clustering quality, significantly improving the biological interpretability and stability of the analysis results.

[0130] The Gaussian kernel function adopts the form of radial basis function, and calculates the nonlinear similarity between samples through the combination of exponential function and squared Euclidean distance. The core principle of the function is based on the probability density function of Gaussian distribution, and the linear inseparable problem in high-dimensional space is transformed into a linear separable problem in kernel space by using distance decay effect. The bandwidth parameter σ kernel controls the decay rate of similarity, and smaller bandwidth value produces more localized similarity measure, and larger bandwidth value produces more global similarity pattern. Compared with the traditional linear similarity measure method, the Gaussian kernel function can better capture the complex nonlinear relationship in high-dimensional molecular data, and provide more accurate similarity information for spectral clustering algorithm, thereby significantly improving the accuracy and robustness of molecular subtype identification.

[0131] The standardized Laplacian matrix solves the singularity problem of the original Laplacian matrix through the standardization of the degree matrix. The construction process of the matrix includes three steps of degree matrix calculation, Laplacian matrix construction and standardization. The standardization of the matrix is realized through the left and right multiplication operation of M -1 / 2 The characteristic decomposition process solves the generalized eigenvalue problem, and the obtained eigenvectors contain important topological information of the graph structure. The eigenvectors corresponding to the smallest eigenvalues reveal the clustering structure of the data. Compared with the traditional K-means clustering method, the spectral clustering algorithm can process clustering structures of any shape, especially suitable for processing complex non-convex clustering patterns in high-dimensional molecular data, and significantly improves the accuracy and scope of disease molecular subtype identification.

[0132] The conversion rate and contribution rate calculation formula quantifies the stability and importance of the classification result through a simple ratio form. The conversion rate reflects the stability of the classification result by counting the proportion of the number of molecules that have changed classification to the total number. When the conversion rate exceeds the preset threshold, the parameter re-optimization mechanism is triggered. The contribution rate evaluates the importance of different expression classifications on the final result by calculating the change degree of the clustering result before and after removing the expression classification. Compared with the static quality control method, the dynamic evaluation mechanism has stronger self-adaptability and error correction function, can timely discover and correct unstable factors in the analysis process, and ensures the reliability and reproducibility of the molecular subtype identification result, providing more reliable molecular diagnostic basis for clinical application.

[0133] In order to better understand and implement the present application, the following provides an embodiment 2 of a specific application scenario of the present application: a technical team uses the method of the present application to identify molecular subtypes of lung cancer patients, in order to provide technical support for precision medical diagnosis. The technical team collected 180 blood samples of lung cancer patients and 60 blood samples of healthy control group, the sample sources were all morning fasting venous blood, the blood volume was 12 milliliters, including 89 cases of adenocarcinoma, 56 cases of squamous cell carcinoma and 35 cases of small cell lung cancer.

[0134] The technical team strictly performed the sample collection and data acquisition process according to step S01. After the blood sample was placed at room temperature for 30 minutes, it was centrifuged at 3000 revolutions per minute for 15 minutes at 4℃, and the serum and plasma components were successfully separated. Illumina HiSeq X Ten platform was used for whole transcriptome sequencing, and the average gene expression profile data of 42 million reads per sample was obtained. Serum proteome was detected by LC-MS / MS technology, and 1847 proteins were identified. A total of 326 metabolites were detected by GC-MS technology for plasma metabolome analysis.

[0135] In the quality control phase of step S02, the technical team detected the hemoglobin concentration range of 118-165 grams per liter, which met the normal range requirements. The degree of hemolysis evaluation showed that the free hemoglobin concentration of all samples was less than 0.15 grams per liter, and there was no hemolysis phenomenon. The storage time verification showed that the time from blood collection to processing was controlled within 1.5 hours, and the frozen storage temperature was stable at -80℃. DNA integrity detection was evaluated by agarose gel electrophoresis, and all samples showed clear high molecular weight bands, with a quality pass rate of 100%.

[0136] In step S03, the technical team used phenol chloroform extraction method for nucleic acid extraction, and the purity ratio of DNA and RNA obtained was between 1.85-2.02. Protein separation used ammonium sulfate fractionation, and after ion exchange chromatography purification, the protein purity reached more than 97%. Gel filtration chromatography used Sephadex G-75 gel medium to effectively separate protein components of different molecular weights.

[0137] In step S04, the technical team established a sequencing quality evaluation function F seq = 0.3·D depth + 0.25·L read + 0.25·Q score + 0.2·C coverage . The initial sequencing parameters were set as: sequencing depth 30 times, read length 150 base pairs, quality score threshold 20, coverage 95%. The initial F seq value was calculated to be 0.68, which was lower than the first threshold value 0.75, so the sequencing depth was adjusted to 50 times. After recalculating, the F seq value increased to 0.82, meeting the quality requirements. The comparison before and after adjusting the sequencing parameters is shown in Figure 4 .

[0138] In step S05, the technical team used ComBat algorithm for batch effect correction, eliminating the systematic differences between the 3 different experimental batches. The results of the intragroup correlation coefficient calculation of technical repeats showed that the correlation coefficient of gene expression data was 0.94, the correlation coefficient of proteome data was 0.91, and the correlation coefficient of metabolome data was 0.89, all meeting the consistency requirements. The DESeq2 method was used for gene expression data normalization, the total protein standardization method was used for proteome data processing, and the internal standardization method was used for metabolome data correction.

[0139] In step S06, the technical team used Z-score standardization method to process the corrected data, and the calculation formula was The missing values were filled by K-Nearest Neighbors imputation algorithm with K value set to 5, successfully filling 2.3% of the missing data. The outliers were detected by the interquartile range method, identifying and processing 0.8% of the outliers. The expression level classification results are shown in Table 1:

[0140] Table 1 Classification statistics of molecular marker expression levels

[0141] Expression classification Number of genes Number of proteins Number of metabolites Total High expression classification 4523 462 82 5067 Medium expression classification 9046 923 163 10132 Low expression classification 4531 462 81 5074

[0142] In step S07, the technical team builds a double-layer game optimization model. The adjustment coefficients in the objective function of the upper-layer model are set to λ = 0.1, μ = 0.05, v = 0.08, and η = 1.5, and the number of clusters n is initially set to 6. The adjustment coefficients in the objective function of the lower-layer model are set to ρ = 0.6, ξ = 0.2, τ = 0.1, and ω = 0.1. After 25 rounds of iterative optimization, the clustering error of the upper-layer model converges to 0.0847, and the feature retention rate of the lower-layer model reaches 0.9234. The convergence process of the double-layer game optimization model is shown in Figure 2 .

[0143] In step S08, the ladder hierarchical sparse coding method is used for feature selection. The technical team sets 9 sparsity levels, decreasing from 0.9 to 0.1. At the sparsity level of 0.9, 2027 features are retained, at the sparsity level of 0.8, 1621 features are further screened, and finally at the sparsity level of 0.3, 486 key features are obtained, with a feature reproducibility rate of more than 85%. The feature selection process of the ladder hierarchical sparse coding is shown in Figure 6 .

[0144] In step S09, the technical team converts the grid optimization problem into a maximum weight matching problem of a bipartite graph and solves it using the Hungarian algorithm. The grid is divided into an 8x8 two-dimensional structure, with a total of 64 grid nodes. The algorithm execution time is 0.24 seconds, and the globally optimal grid division scheme is obtained, with a matching efficiency of 97.8%.

[0145] In step S10, the technical team constructs a similarity matrix using the Gaussian kernel function The bandwidth parameter σ of the kernel function kernel is set to 0.045, which is calculated based on the 486-dimensional feature space. The similarity matrix is a 240x240 symmetric positive definite matrix, with diagonal elements equal to 1 and non-diagonal elements ranging from 0.023 to 0.876. The sample similarity distribution is shown in Figure 5 .

[0146] In step S11, the spectral clustering algorithm is used for unsupervised clustering analysis. The calculation of the degree matrix M shows that the sample degree value distribution ranges from 156.7 to 203.4. The normalized Laplacian matrix L normThe eigenvalue decomposition of the feature matrix obtains 240 eigenvalues, and the first 6 eigenvalues in ascending order are 0.0012, 0.0089, 0.0156, 0.0234, 0.0789, and 0.1245. The eigenvalue gap analysis shows that the difference between the fourth and fifth eigenvalues is the largest, and the optimal cluster number h = 4 is determined. The eigenvalue distribution and cluster number determination process are shown in Figure 8 .

[0147] The cluster result analysis is shown in Table 2:

[0148] Table 2: Lung cancer molecular subtype cluster result statistics

[0149] Subtype number Sample number Adenocarcinoma patients Squamous carcinoma patients Small cell lung cancer patients Healthy controls Subtype 1 67 45 18 4 0 Subtype 2 58 32 21 5 0 Subtype 3 55 12 17 26 0 Subtype 4 60 0 0 0 60

[0150] In step S12, the technical team verifies the biological function of the clustering result. Gene enrichment analysis shows that subtype 1 is significantly enriched in cell proliferation-related pathways, with a P value of 1.2 x 10 -8 ; subtype 2 is enriched in DNA repair pathways, with a P value of 3.4 x 10 -7 ; and subtype 3 is enriched in cell apoptosis pathways, with a P value of 2.1 x 10 -9 . Protein-protein interaction network analysis identifies 23 key regulatory proteins, including important cancer genes such as TP53, EGFR, and KRAS. Metabolic pathway analysis shows that different subtypes have significant differences in glycolysis, fatty acid metabolism, and amino acid metabolism.

[0151] In step S13, the technical team establishes a subtype stability evaluation mechanism. The transition rate calculation result shows that the transition rate of high expression classification is 8.3%, the transition rate of medium expression classification is 12.1%, and the transition rate of low expression classification is 9.7%, all of which are lower than the preset threshold of 15%, indicating that the classification result is stable and reliable. Contribution rate analysis shows that the contribution rate of high expression classification to the clustering result is 34.2%, the contribution rate of medium expression classification is 41.6%, and the contribution rate of low expression classification is 24.2%. The clustering result stability evaluation is shown in Figure 7 .

[0152] In step S14, the technical team constructs a subtype-specific molecular marker spectrum. The specific markers of subtype 1 include 67 genes, 12 proteins, and 8 metabolites, among which the EGFR gene expression fold change is 3.7 times, and the P value is 2.3 x 10 -12The specific markers of subtype 2 include 54 genes, 9 proteins and 6 metabolites, and the expression fold change of BRCA1 gene is 4.2 times. The specific markers of subtype 3 include 73 genes, 15 proteins and 11 metabolites, and the expression fold change of BAX protein is 5.1 times. The diagnostic performance evaluation shows that the area under the receiver operating characteristic curve of the combination of molecular markers reaches 0.923, which has excellent diagnostic performance. The number of specific markers of each molecular subtype is compared as shown in Table 3. Figure 3

[0153] The final identification result is shown in Table 3:

[0154] Table 3: Number of key molecular markers of each subtype

[0155] Subtype number Number of key genes Number of key proteins Number of key metabolites Diagnostic accuracy Subtype 1 67 12 8 94.3％ Subtype 2 54 9 6 91.7％ Subtype 3 73 15 11 96.4％ Healthy controls 41 7 4 100.0％

[0156] The present application brings significant technical progress compared to traditional molecular subtype identification methods. Traditional methods usually use single omics data for analysis, which is easily affected by data noise and batch effects. The present application integrates multi-omics data for analysis and constructs a more comprehensive and stable molecular feature map. The application of the double-layer game optimization model solves the problem that traditional single-objective optimization methods are prone to local optimization when dealing with complex high-dimensional data. Through the coordinated optimization of upper and lower layers, the best balance between clustering accuracy and feature retention rate is achieved. The ladder grading sparse coding method can better capture the hierarchical relationship and dependency between features compared to the traditional one-time feature selection strategy, avoiding the loss of important biological information. The establishment of the dynamic stability evaluation mechanism enables the entire analysis process to have adaptive adjustment capability, which can automatically trigger parameter re-optimization when detecting unstable classification results, ensuring the reliability and reproducibility of the analysis results. The application of the spectral clustering algorithm combined with the Gaussian kernel function enables the method to handle clustering structures of any shape, especially suitable for non-convex clustering patterns commonly found in molecular data, significantly improving the accuracy of complex molecular subtype identification.

[0157] It should be noted that the variables involved in the present application are explained in detail as shown in Table 4.

[0158] Table 4: Variable explanation table

[0159]

[0160]

[0161] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited thereto, and any person skilled in the art can easily think of changes or replacements within the technical scope disclosed by the present application, which should be covered within the protection scope of the present application.​

Claims

1. A method of unsupervised clustering analysis of disease marker molecular subtypes, characterized in that, The application relates to a method for constructing a disease subtype molecular marker spectrum. Blood samples of patients with diseases and blood samples of healthy controls are collected, and original molecular marker data containing gene expression profile data, proteomic data and metabolomic data are obtained through high-throughput sequencing technology; Quality control detection and nucleic acid extraction and protein separation are carried out on the blood samples; a sequencing quality evaluation function is established to optimize sequencing parameters; Batch effect correction and preprocessing are carried out on the original molecular marker data to obtain a standardized molecular marker data matrix classified into high expression classification, medium expression classification and low expression classification according to expression levels; a double-layer game optimization model is established to optimize grid division parameters; the upper model takes minimizing clustering error as an objective, the lower model takes maximizing feature retention rate as an objective, and the two objective functions are associated through a coupling term; based on the standardized molecular marker data matrix, a ladder hierarchical sparse coding method is adopted for feature selection; The grid optimization problem is constructed into a maximum weight matching problem of a bipartite graph, and the optimal grid division parameters are solved through the Hungarian algorithm; a similarity matrix is constructed, and unsupervised clustering analysis is carried out through a spectral clustering algorithm; The clustering results are verified in biological functions; a subtype stability evaluation mechanism is established; and a subtype-specific molecular marker spectrum is established.

2. The disease marker molecule subtype unsupervised clustering analysis method according to claim 1, characterized by, The high-throughput sequencing technology refers to a technology for simultaneously detecting tens of thousands to millions of molecular markers in a single experiment through parallel sequencing of a large number of DNA or RNA fragments.

3. The disease marker molecule subtype unsupervised clustering analysis method according to claim 2, characterized by, The quality control detection specifically includes hemoglobin concentration determination, hemolysis degree evaluation and storage time verification, and blood samples meeting the quality standards are screened out for subsequent analysis.

4. The disease marker molecule subtype unsupervised clustering analysis method according to claim 3, characterized by, The nucleic acid extraction and protein separation specifically adopts column chromatography and gel electrophoresis to purify the extract to obtain a high-purity molecular marker extract.

5. The disease marker molecule subtype unsupervised clustering analysis method according to claim 4, characterized by, The column chromatography refers to a chromatographic technology for separating a mixture by utilizing the interaction difference between a stationary phase and a mobile phase.

6. The disease marker molecule subtype unsupervised clustering analysis method according to claim 5, characterized by, The gel electrophoresis refers to a technical method for separating by utilizing the size and charge difference under the action of an electric field.

7. The disease marker molecule subtype unsupervised clustering analysis method according to claim 6, characterized by, The sequencing quality evaluation function specifically comprehensively evaluates the sequencing data quality; when the evaluation function value is less than a first threshold value, the sequencing depth parameter is adjusted; when the evaluation function value is less than a second threshold value, the read length parameter is adjusted; and when the evaluation function value is less than a third threshold value, the quality threshold value parameter is improved.

8. The disease marker molecule subtype unsupervised clustering analysis method according to claim 7, characterized by, The batch effect correction refers to a data processing method for eliminating systematic differences between different experimental batches; the systematic error between different batches is eliminated through a standardization method to obtain corrected molecular marker data.

9. The disease marker molecule subtype unsupervised clustering analysis method according to claim 8, characterized by, The high expression classification refers to a molecular marker category with an expression level above the upper quartile of data distribution; the medium expression classification refers to a molecular marker category with an expression level near the median of data distribution; and the low expression classification refers to a molecular marker category with an expression level below the lower quartile of data distribution.

10. The disease marker molecule subtype unsupervised clustering analysis method according to claim 9, characterized by, The ladder hierarchical sparse coding method refers to a coding method for gradually screening important features from coarse to fine by setting multiple sparsity levels; key features are screened level by level by setting different sparsity threshold values to obtain a feature subset after dimension reduction.

Images