Data classification method and system for single cell sequencing

Abstract

The invention relates to the technical field of bioinformatics, in particular to a data classification method and system for single cell sequencing. The method comprises the following steps: acquiring an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; constructing a cell topological graph according to the preprocessed gene expression matrix to obtain the cell topological graph; performing cell division rate calculation according to the cell topological graph and the preprocessed gene expression matrix, and performing RNA rate calculation to obtain RNA rate data; adding node rates to the cell topological graph according to the cell division rate data and the RNA rate data to obtain an enhanced cell topological graph. According to the method, the dynamic graph is constructed and the change rule of the dynamic graph is analyzed, so that the dynamic change of the cell state is better captured, the sensitivity to parameters is reduced, a complex cell population structure is processed, and a quantitative cell heterogeneity evaluation index is provided.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a data classification method and system for single-cell sequencing. Background Art

[0002] Early single-cell sequencing data analysis mainly relied on traditional clustering methods, such as k-means and hierarchical clustering. These methods group cells based on the similarity of gene expression, but they usually assume that cell populations are static and ignore the dynamic changes in cell states. With the advancement of single-cell sequencing technology and the increase in data volume, new analysis methods have gradually emerged, such as graph-based clustering methods, dimensionality reduction and visualization-based methods, and trajectory inference methods. However, existing methods are insufficient in analyzing cell heterogeneity.

[0003] Insufficient analysis of cell heterogeneity:

[0004] Difficult to capture the dynamic changes of cell states: Traditional clustering methods divide cells into fixed groups, ignoring the continuity and dynamic changes of cell states.

[0005] Sensitive to clustering algorithm parameters: Clustering results are usually sensitive to the choice of algorithm parameters, such as the k value in k-means. Different parameter choices will lead to different clustering results, affecting the understanding of cell heterogeneity.

[0006] Difficulty in handling complex cell population structures: Some cell populations have complex hierarchical structures or continuous state changes, and traditional clustering methods are difficult to handle these situations effectively.

[0007] Lack of quantitative assessment of cell heterogeneity: Many methods can only qualitatively describe cell heterogeneity and lack quantitative assessment indicators for cell heterogeneity. Summary of the invention

[0008] Based on this, it is necessary to provide a single-cell sequencing data classification method and system to solve at least one of the above technical problems.

[0009] To achieve the above purpose, a single cell sequencing data classification method comprises the following steps:

[0010] Step S1: obtaining an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; constructing a cell topology map according to the preprocessed gene expression matrix to obtain a cell topology map;

[0011] Step S2: Calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights according to the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference according to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights to obtain a dynamic trajectory map;

[0012] Step S3: Perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix;

[0013] Step S4: performing clustering probability distribution calculation according to the cell membership matrix and the preprocessed gene expression matrix, and performing gene sorting to obtain gene entropy sorting data; performing module information entropy calculation according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; performing feature gene screening according to the gene entropy sorting data and the module information entropy sorting data to obtain a feature gene list;

[0014] Step S5: Annotate the cell state according to the characteristic gene list to obtain the cell state annotation result to realize the data classification task of single-cell sequencing.

[0015] The present invention effectively removes noise data, highlights highly variable genes, and constructs a topological map based on local clustering results by preprocessing and constructing a cell topological map, more accurately reflects the connection relationship between cells and the local structure in high-dimensional space, and provides a more reliable and more sophisticated basis for subsequent analysis. By calculating and integrating cell division rate and RNA rate information, an enhanced cell topological map is constructed, and the dynamic trajectory map is inferred using integrated weights and optimal path algorithms, effectively capturing the dynamic changes of cells and the dynamic trends of gene expression, more accurately depicting the cell development trajectory, and providing dynamic evolution information for cell state annotation. By performing fuzzy cluster analysis on the dynamic trajectory map and constructing a cell membership matrix, cell subpopulations are effectively identified, and the degree to which each cell belongs to different clusters is quantified, providing more refined and accurate classification information for subsequent feature gene screening and cell state annotation. By calculating gene information entropy and module information entropy, and combining the cell membership matrix for weighting and sorting, genes and functional modules specifically expressed in different cell subpopulations are effectively identified, thereby screening out the feature genes that contribute the most to cell state annotation, and improving the accuracy and biological significance of cell state annotation. By using a preset cell marker gene database and screened characteristic genes, training and optimizing a multi-label classification model, and predicting and annotating cell states, the automatic classification of single-cell sequencing data is achieved, and a more comprehensive and accurate cell state annotation result is provided, which facilitates researchers to conduct in-depth analysis and interpretation of single-cell data. Therefore, the present invention provides a data classification method for cell sequencing. By constructing a dynamic graph and analyzing its changing rules, the method can better capture the dynamic changes of cell states, reduce sensitivity to parameters, process complex cell population structures, and provide quantitative cell heterogeneity evaluation indicators, thereby overcoming the shortcomings of existing methods for cell heterogeneity analysis and deepening the understanding of cell heterogeneity.

[0016] Preferably, step S1 comprises the following steps:

[0017] Step S11: obtaining an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix;

[0018] Step S12: performing high-dimensional space distance calculation on the distances between cells in the preprocessed gene expression matrix to obtain high-dimensional space distance data;

[0019] Step S13: selecting a filter according to the high-dimensional space distance data to obtain a filter function; dividing the value range of the filter function into multiple overlapping intervals according to the high-dimensional space distance data to obtain interval division data;

[0020] Step S14: performing local clustering on the preprocessed gene expression matrix according to the interval division data to obtain a local clustering result;

[0021] Step S15: constructing a cell topology map according to the local clustering results to obtain a cell topology map.

[0022] The present invention can remove noise and low-quality data, such as low-expression genes, dead cells or background RNA, by obtaining the original gene expression matrix and performing preprocessing, thereby improving the accuracy and reliability of subsequent analysis. The normalization step eliminates the difference in sequencing depth between different cells, the logarithmic transformation can stabilize the data variance, and the selection of highly variable genes focuses on the genes that contribute the most to cell heterogeneity, making the constructed cell topology more biologically meaningful. By calculating the distance between cells in high-dimensional space, the transcriptome similarity between cells can be quantified. This provides a basis for the subsequent construction of cell topology, so that the topology can accurately reflect the relationship between cells. Selecting a suitable filter function and dividing its value range into multiple overlapping intervals can effectively capture the local neighboring relationship between cells. The design of overlapping intervals allows cells to belong to multiple local neighborhoods at the same time, thereby better reflecting the continuous distribution of cells in high-dimensional space, and avoiding the information loss caused by hard division. Local clustering is performed on each overlapping interval, and the local structure between cells at different scales can be identified. This helps to build a more sophisticated cell topology, which can reflect the multi-level relationship between cells, rather than just based on global similarity. Constructing a cell topology map based on local clustering results can more accurately reflect the connection relationship between cells and retain the local structural information of cells in high-dimensional space. Compared with directly constructing a topology map based on the distance matrix, this method can better handle the noise and complex structure in high-dimensional data, making the constructed topology map more robust and better reflecting the real cell relationship.

[0023] Preferably, step S2 comprises the following steps:

[0024] Step S21: Calculate the cell division rate according to the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculate the RNA rate of each cell using the pre-processed gene expression matrix to obtain RNA rate data;

[0025] Step S22: adding node attributes to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map;

[0026] Step S23: Initializing the node weights according to the enhanced cell topology map to obtain initial node weights;

[0027] Step S24: adjusting the node weights of the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain the division rate adjustment weights; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain the RNA rate adjustment edge weights;

[0028] Step S25: integrating the node weights and edge weights of the division rate adjustment weights and the RNA rate adjustment edge weights according to the enhanced cell topology map to obtain an integrated weight;

[0029] Step S26: applying an optimal path algorithm according to the enhanced cell topology map and the integrated weights to obtain optimal path data;

[0030] Step S27: constructing a dynamic trajectory map for the enhanced cell topology map according to the optimal path data to obtain a dynamic trajectory map.

[0031] The present invention calculates the cell division rate and RNA rate respectively by combining the cell cycle related gene expression and the mRNA information before and after splicing, and can more accurately capture the dynamic changes of cells. The cell division rate reflects the active degree of cell proliferation, while the RNA rate reflects the dynamic change trend of gene expression. These two indicators provide important biological information for the subsequent dynamic trajectory inference. The cell division rate and RNA rate are added as node attributes to the cell topology map, and an enhanced cell topology map is constructed, so that the topology map not only contains the connection information between cells, but also contains the dynamic information of each cell itself. This provides a basis for the subsequent node weight and edge weight adjustment, so that the dynamic trajectory inference can better reflect the real process of cell development. Initialize the node weight, and provide a starting point for the subsequent weight adjustment based on the cell division rate and RNA rate. The initialization weight is 1, which ensures that all cells have the same weight in the initial state, avoiding the influence of the initial value on the subsequent calculation results. The node weight is adjusted based on the cell division rate, so that cells with high division rate obtain higher weights in trajectory inference, highlighting the importance of these cells in the cell development process. At the same time, the edge weights are adjusted based on the RNA rate, so that the edges with large RNA rate changes obtain higher weights, thereby better capturing the impact of dynamic changes in gene expression on cell state transitions. The division rate adjustment weights and RNA rate adjustment edge weights are integrated to comprehensively consider the impact of cell division and dynamic changes in gene expression on cell state transitions, so that the integrated weights can better reflect the driving force of cell development. The optimal path algorithm is applied to the enhanced cell topology map to identify the most likely paths connecting different cell states, which represent the trajectory of cell development. Using the integrated weights as the basis for path calculation, the optimal path can reflect the combined impact of cell division and dynamic changes in gene expression. A dynamic trajectory map is constructed based on the optimal path data, which converts the static connection relationship between cells into a dynamic evolution trajectory, more intuitively showing the cell state transition process, and providing an important basis for subsequent cell state annotation and analysis.

[0032] Preferably, step S21 includes the following steps:

[0033] Step S211: extracting the cell cycle-related gene expression matrix from the preprocessed gene expression matrix to obtain a cycle gene expression matrix;

[0034] Step S212: performing single cell cycle scoring according to the cycle gene expression matrix to obtain a single cell cycle score;

[0035] Step S213: dividing the cell cycle stages according to the single cell cycle scores to obtain the cell cycle stages;

[0036] Step S214: Calculate the cell division rate according to the cell topology and the cell cycle stage to obtain the cell division rate;

[0037] Step S215: constructing a pre- and post-splicing mRNA matrix for the pre-processed gene expression matrix to obtain a pre- and post-splicing mRNA matrix;

[0038] Step S216: Infer the RNA rate based on the mRNA matrix before and after splicing to obtain the RNA rate matrix; perform rate smoothing and noise reduction on the RNA rate matrix to obtain RNA rate data.

[0039] The present invention extracts the cell cycle-related gene expression matrix, focuses on genes closely related to the cell cycle process, removes the interference of other gene expression information, makes the cell cycle score more accurate, and can better reflect the proliferation state of the cell. By calculating the single cell cycle score, the position of each cell in the cell cycle can be quantified, providing a basis for the subsequent cell cycle stage division. This can more accurately reflect the cell cycle state than relying solely on a few marker genes, because it takes into account the comprehensive expression pattern of multiple genes. Dividing cells into different cell cycle stages (G1 / S / G2 / M) can more finely describe the proliferation state of the cell. This is crucial for understanding cell heterogeneity and developmental trajectories, because cells in different cycle stages have different transcriptome characteristics and biological functions. The cell division rate is calculated in combination with the cell topology map and the cell cycle stage information, taking into account the relationship between the cell and its neighbor cells, and can more accurately reflect the local proliferation activity of the cell. For example, even if a cell has a high cycle score, if its neighbor cells are all in the stationary phase, then its actual division rate may not be high. Constructing the mRNA matrix before and after splicing provides key data for RNA rate inference. Pre-splicing mRNA represents newly synthesized transcripts, while post-splicing mRNA represents mature transcripts. By comparing the ratio of the two, the dynamic trend of gene expression can be inferred. Inferring RNA rate based on the pre- and post-splicing mRNA matrix can capture the dynamic changes of gene expression, and remove the impact of technical errors through smoothing and noise reduction, making RNA rate data more reliable and more able to reflect the true dynamics of gene expression. This provides important information for studying cell state transitions and developmental trajectories.

[0040] Preferably, step S24 includes the following steps:

[0041] Step S241: extracting the average node division rate of the enhanced cell topology map to obtain the average node division rate; setting the minimum division rate according to the average node division rate to obtain the minimum division rate;

[0042] Step S242: adjusting the initial node weights according to the average node splitting rate and the minimum node splitting rate to obtain a splitting rate adjustment weight;

[0043] Step S243: extracting the node average RNA rate vector of the enhanced cell topology map to obtain the node average RNA rate vector; performing edge RNA rate difference calculation on the enhanced cell topology map according to the node average RNA rate vector to obtain edge RNA rate difference data;

[0044] Step S244: setting the minimum RNA rate difference according to the edge RNA rate difference data to obtain the minimum RNA rate difference data;

[0045] Step S245: calculating the directional weight adjustment factor of the edge of the enhanced cell topology map according to the edge RNA rate difference data to obtain the directional weight adjustment factor;

[0046] Step S246: adjusting the edge weight of the edge RNA rate difference data according to the minimum RNA rate difference data and the directional weight adjustment factor to obtain the RNA rate adjusted edge weight.

[0047] The present invention provides a benchmark for subsequent node weight adjustment by extracting the average node division rate and setting the minimum division rate. Using the average division rate can avoid the influence of individual cell division rate abnormal values ​​on the results, while setting the minimum division rate (especially the non-zero minimum value) can avoid the situation of division by zero and ensure the proportional relationship of weight adjustment. Adjusting the node weight according to the average node division rate and the minimum division rate highlights the importance of cells with high division rates in dynamic trajectory inference. Cells with high division rates are usually in an active proliferation state and may represent key nodes or transition states of cell development, so giving them higher weights helps to more accurately identify cell development trajectories. Extracting the node average RNA rate vector and calculating the edge RNA rate difference can quantify the degree of difference in the dynamic changes of gene expression between cells. This provides a basis for subsequent edge weight adjustment, so that the edge weight can reflect the amplitude of gene expression changes during cell state conversion. Setting the minimum RNA rate difference is similar to setting the minimum division rate in step S241, which can avoid the situation of division by zero and ensure the proportional relationship of edge weight adjustment. This makes the weight adjustment more robust and less susceptible to extreme values. The directional weight adjustment factor for edges is calculated, taking into account the relative size of RNA rate differences, further highlighting the importance of edges with large RNA rate changes in dynamic trajectory inference. This helps identify cell state transition events with dramatic gene expression changes. The edge weights are adjusted based on the minimum RNA rate difference and the directional weight adjustment factor, so that edges with large RNA rate changes receive higher weights. This highlights the driving effect of dynamic changes in gene expression on cell state transitions, making the constructed dynamic trajectory graph more reflective of the real process of cell development.

[0048] Preferably, step S3 comprises the following steps:

[0049] Step S31: performing feature matrix conversion on the dynamic trajectory graph to obtain a trajectory feature matrix;

[0050] Step S32: determining the number of clusters for the trajectory feature matrix to obtain cluster number data;

[0051] Step S33: Perform fuzzy C-means clustering according to the trajectory feature matrix and the cluster quantity data to obtain a fuzzy clustering result;

[0052] Step S34: constructing a membership matrix according to the fuzzy clustering results to obtain a cell membership matrix.

[0053] The present invention converts the dynamic trajectory graph into a trajectory feature matrix and converts the graph structure information into a numerical matrix, which facilitates subsequent cluster analysis. Feature decomposition extracts the main features of the dynamic trajectory graph, reduces the dimension of the data, and retains key information, making the clustering results more robust and easier to interpret. The silhouette coefficient method is used to determine the number of clusters, avoiding the deviation caused by the subjective selection of the number of clusters, making the clustering results more objective and reasonable. The silhouette coefficient can effectively evaluate the quality of the clustering results. Selecting the K value with the highest silhouette coefficient can find the optimal number of clusters, so as to better identify cell subpopulations. Clustering the trajectory feature matrix using the fuzzy C-means clustering algorithm can effectively handle the transitional states that may exist between cells in different states. Fuzzy clustering allows a cell to belong to multiple clusters at the same time and assign different memberships, which can better reflect the continuity and dynamics of the cell state than hard clustering. The cell membership matrix is ​​constructed to quantify the degree to which each cell belongs to different clusters. This provides important information for subsequent feature gene screening and cell state annotation, for example, the cell type or state to which it belongs can be judged according to the membership of the cell.

[0054] Preferably, step S4 comprises the following steps:

[0055] Step S41: clustering the preprocessed gene expression matrix according to the cell membership matrix to obtain a gene expression grouping list;

[0056] Step S42: performing clustering probability distribution calculation according to the gene expression grouping list, and performing gene sorting to obtain gene entropy sorting data;

[0057] Step S43: selecting a co-expression measurement method according to the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters;

[0058] Step S44: Calculate the inter-gene co-expression values ​​of the pre-processed gene expression matrix using the co-expression measurement method and parameters to obtain a gene co-expression matrix;

[0059] Step S45: constructing a co-expression network according to the gene co-expression matrix to obtain a gene co-expression network; identifying network modules of the gene co-expression network to obtain network modules;

[0060] Step S46: Calculate the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain the module information entropy ranking data;

[0061] Step S47: Screen characteristic genes according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list.

[0062] The present invention can analyze the gene expression pattern for each cell subgroup separately by grouping the gene expression matrix by clustering, thereby identifying the genes specifically expressed in different subgroups. This helps to understand the differences between different cell types or states, and provides a basis for subsequent characteristic gene screening. By calculating the gene clustering probability distribution and gene information entropy, the specificity of gene expression in different cell subgroups can be quantified. The higher the information entropy, the greater the difference in gene expression in different subgroups, and the more likely it is a characteristic gene. The weighted gene information entropy further considers the cell membership, making the selection of characteristic genes more accurate. Selecting appropriate co-expression measurement methods and parameters can accurately evaluate the co-expression relationship between genes. Pearson correlation coefficient is a commonly used co-expression measurement method that can capture the linear relationship between gene expression. Appropriate parameter settings can effectively filter out weakly correlated or falsely correlated gene pairs. Calculate the co-expression value between genes, construct a gene co-expression matrix, and provide a data basis for the subsequent construction of a gene co-expression network. The gene co-expression matrix reflects the similarity of expression patterns between genes and can be used to identify functionally related gene modules. Constructing a gene co-expression network and identifying network modules can reveal the interaction relationship and functional modules between genes. Modules in co-expression networks are usually composed of functionally related genes, which can reflect biological pathways or regulatory mechanisms within cells. Calculating the information entropy of network modules can evaluate the specificity of modules in different cell subsets. The higher the information entropy, the greater the expression difference of the module in different subsets, and it is more likely to be a functional module related to a specific cell type or state. Combining gene entropy sorting data and module information entropy sorting data for feature gene screening can more effectively identify feature genes related to cell state annotations. This method comprehensively considers the individual specificity and module specificity of gene expression, making the screened feature genes more representative and biologically meaningful.

[0063] Preferably, step S42 includes the following steps:

[0064] Step S421: determining a kernel function according to the gene expression grouping list to obtain a kernel function; determining a bandwidth according to the gene expression grouping list and the kernel function to obtain a bandwidth parameter;

[0065] Step S422: performing kernel density estimation on the gene expression grouping list according to the kernel function and the bandwidth parameter to obtain the gene clustering probability density; performing probability density normalization on the gene clustering probability density to obtain the gene clustering probability distribution data;

[0066] Step S423: Calculate gene information entropy according to gene clustering probability distribution data to obtain gene information entropy;

[0067] Step S424: adjusting the information entropy weight of the gene information entropy according to the cell membership matrix to obtain weighted gene information entropy;

[0068] Step S425: performing gene entropy sorting according to the weighted gene information entropy to obtain gene entropy sorting data.

[0069] The present invention can effectively perform kernel density estimation by selecting a Gaussian kernel function and using the Silverman's rule of thumb method to determine the bandwidth. The Gaussian kernel function has good mathematical properties, and the Silverman's rule of thumb method can adaptively select the bandwidth according to the distribution characteristics of the data, so that the result of the kernel density estimation is more accurate. The kernel density estimation method is used to calculate the probability density of gene clustering, which can more accurately describe the distribution of gene expression in different cell subpopulations. Compared with the simple histogram method, the kernel density estimation can better handle the noise and fluctuation in the data and provide a smoother probability density curve. The probability density normalization ensures the validity of the probability distribution. The gene information entropy is calculated to quantify the degree of dispersion of gene expression in different cell subpopulations. The higher the information entropy, the greater the difference in gene expression in different subpopulations, and it is more likely to be a characteristic gene. Using the cell membership matrix to weight the gene information entropy can more accurately reflect the importance of genes in different cell subpopulations. Considering the cell membership can avoid the information loss caused by fuzzy clustering and highlight those genes that are highly expressed in specific subpopulations and have high information entropy. By sorting genes according to weighted gene information entropy, genes with high information entropy and high expression in specific cell subpopulations can be ranked first. This helps in subsequent feature gene screening and gives priority to genes that contribute most to cell state annotation.

[0070] Preferably, step S5 comprises the following steps:

[0071] Step S51: preparing training data according to the characteristic gene list and the preset cell marker gene database to obtain a training data set;

[0072] Step S52: training a multi-label classification model using the training data set to obtain a multi-label classification model;

[0073] Step S53: performing model evaluation and optimization on the multi-label classification model to obtain an optimized multi-label classification model;

[0074] Step S54: According to the preprocessed gene expression matrix and the characteristic gene list, the optimized multi-label classification model is used to predict the cell state to obtain the cell state prediction result;

[0075] Step S55: Generate a cell state annotation result according to the cell state prediction result to obtain a cell state annotation result.

[0076] The present invention prepares a training data set by using a characteristic gene list and a preset cell marker gene database, and provides high-quality annotation data for the training of a multi-label classification model. Reliable cell type marker gene information can be obtained by using the CellMarker database, and a training data set rich in information and biological significance can be constructed in combination with the screened characteristic genes. The multi-label classification model is trained using the training data set so that the model can learn the association between the characteristic gene expression pattern and the cell type. The RandomForestClassifier algorithm has good generalization ability and robustness, and is suitable for processing high-dimensional data and multi-label classification problems. Through model evaluation and optimization, the prediction performance of the multi-label classification model can be improved. The generalization ability of the model is evaluated using cross validation and Hamming loss, and the best model parameters are found using grid search, which can effectively avoid overfitting and improve the prediction accuracy of the model. The optimized multi-label classification model is used to predict the cell state of the preprocessed gene expression matrix, and the cell type or state to which each cell belongs can be predicted according to the expression pattern of the characteristic gene. This provides a quantitative prediction result for subsequent cell state annotation. The cell state annotation result is generated according to the cell state prediction result, and the prediction result is converted into an easy-to-understand cell type or state label. This allows researchers to conveniently classify and interpret single-cell data and gain insight into cellular heterogeneity and functional diversity.

[0077] Preferably, the present invention further provides a single-cell sequencing data classification system, which is used to perform the single-cell sequencing data classification method as described above, and the single-cell sequencing data classification system comprises:

[0078] A manifold embedding module is used to obtain an original gene expression matrix; perform gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; and construct a cell topology map according to the preprocessed gene expression matrix to obtain a cell topology map;

[0079] The dynamic trajectory inference module is used to calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights according to the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference according to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights to obtain a dynamic trajectory map;

[0080] The fuzzy clustering module is used to perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix;

[0081] The characteristic gene screening module is used to calculate the clustering probability distribution according to the cell membership matrix and the preprocessed gene expression matrix, and perform gene sorting to obtain gene entropy sorting data; calculate the module information entropy according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; perform characteristic gene screening according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list;

[0082] The cell state annotation module is used to annotate the cell state according to the characteristic gene list and obtain the cell state annotation result to realize the data classification task of single-cell sequencing. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 A schematic diagram of the steps of a data classification method for single-cell sequencing;

[0084] Figure 2 for Figure 1 Detailed implementation steps of step S2 in the flowchart;

[0085] Figure 3 for Figure 1 Detailed implementation steps of step S4 in FIG.

[0086] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings in conjunction with the embodiments. DETAILED DESCRIPTION

[0087] The technical method of the present invention is described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by technicians in this field without creative work are within the scope of protection of the present invention.

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

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

[0090] To achieve this, please refer to Figures 1 to 3 , a data classification method for single-cell sequencing, comprising the following steps:

[0091] Step S1: obtaining an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; constructing a cell topology map according to the preprocessed gene expression matrix to obtain a cell topology map;

[0092] Step S2: Calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights according to the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference according to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights to obtain a dynamic trajectory map;

[0093] Step S3: Perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix;

[0094] Step S4: performing clustering probability distribution calculation according to the cell membership matrix and the preprocessed gene expression matrix, and performing gene sorting to obtain gene entropy sorting data; performing module information entropy calculation according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; performing feature gene screening according to the gene entropy sorting data and the module information entropy sorting data to obtain a feature gene list;

[0095] Step S5: Annotate the cell state according to the characteristic gene list to obtain the cell state annotation result to realize the data classification task of single-cell sequencing.

[0096] In the embodiment of the present invention, reference Figure 1 FIG. 1 is a schematic diagram of the steps of the data classification method for single-cell sequencing of the present invention. In this example, the data classification method for single-cell sequencing includes the following steps:

[0097] Step S1: obtaining an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; constructing a cell topology map according to the preprocessed gene expression matrix to obtain a cell topology map;

[0098] In an embodiment of the present invention, after reading the original gene expression matrix, Scanpy is used for preprocessing, including filtering low-quality cells and low-expression genes, normalizing cell expression, logarithmic transformation, and high-variable gene selection. Then, based on the preprocessed gene expression matrix, the Euclidean distance between cells is calculated, the Gaussian kernel function is selected as a filter, and its range is divided into multiple overlapping intervals. The Louvain algorithm is used for local clustering in each interval, and finally all local clustering results are integrated to construct a cell topology map, and the weight of the edge is the average value of the Gaussian kernel function value of all intervals containing the edge.

[0099] Step S2: Calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights according to the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference according to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights to obtain a dynamic trajectory map;

[0100] In an embodiment of the present invention, based on a cell topology map and a pre-processed gene expression matrix, a cell cycle-related gene expression matrix is ​​extracted, a cell cycle score is calculated, and cells are divided into different cycle stages. The cell division rate is calculated according to the cell cycle stage and the topology map. At the same time, velocity is used to calculate the RNA rate, and smooth noise reduction is performed. The cell division rate and the RNA rate are added to the node attributes of the cell topology map to obtain an enhanced cell topology map. The node weight is initialized to 1, and then the node weight is adjusted according to the cell division rate and the minimum division rate. The edge RNA rate difference is calculated according to the node average RNA rate vector, and the edge weight is adjusted according to the minimum RNA rate difference and the directional weight adjustment factor. Finally, the node weight and the edge weight are integrated, the Dijkstra algorithm is used to calculate the optimal path, and a dynamic trajectory map is constructed according to the optimal path.

[0101] Step S3: Perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix;

[0102] In the embodiment of the present invention, the dynamic trajectory graph is converted into an adjacency matrix, eigendecomposition is performed, and the first K eigenvectors (the K value is determined according to the variance explanation rate) are selected to form a trajectory feature matrix. The silhouette coefficient method is used to determine the optimal number of clusters K, and then the fuzzy C-means clustering algorithm is used to cluster the trajectory feature matrix, and the fuzzy coefficient m is set to 2. Finally, a cell membership matrix is ​​constructed based on the clustering results, and the matrix elements represent the membership of the cell to each cluster.

[0103] Step S4: performing clustering probability distribution calculation according to the cell membership matrix and the preprocessed gene expression matrix, and performing gene sorting to obtain gene entropy sorting data; performing module information entropy calculation according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; performing feature gene screening according to the gene entropy sorting data and the module information entropy sorting data to obtain a feature gene list;

[0104] In an embodiment of the present invention, the gene expression data in the pre-processed gene expression matrix are grouped by clustering according to the cell membership matrix. For each grouping, the bandwidth determined by the Gaussian kernel function and the Silverman's rule of thumb method is used to perform kernel density estimation, calculate the gene clustering probability distribution, and calculate the gene information entropy. Then, the gene information entropy is weighted according to the cell membership matrix. Genes are sorted according to the weighted gene information entropy. At the same time, the co-expression value between genes is calculated using the Pearson correlation coefficient, a gene co-expression network is constructed, and the network module is identified using the Louvain algorithm. The module information entropy of each module is calculated, and finally the characteristic genes are screened according to the gene entropy sorting data and the module information entropy sorting data.

[0105] Step S5: annotating the cell state according to the characteristic gene list to obtain the cell state annotation result to realize the data classification task of single cell sequencing;

[0106] In an embodiment of the present invention, the CellMarker database is used as a preset cell marker gene database to prepare a training data set. The feature vector of the data set consists of the expression values ​​of the feature genes in different cell types, and the label is a binary vector representing the cell type. The multi-label classification model is trained using RandomForestClassifier, and the parameters are set to n_estimators=100, max_depth=10, and random_state=0. The model performance is evaluated using 5-fold cross validation and Hamming loss, and the model parameters are optimized using grid search. The optimized model is used to predict the cell state of the preprocessed gene expression matrix, and finally the cell state annotation results are generated based on the prediction results.

[0107] Preferably, step S1 comprises the following steps:

[0108] Step S11: obtaining an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix;

[0109] Step S12: performing high-dimensional space distance calculation on the distances between cells in the preprocessed gene expression matrix to obtain high-dimensional space distance data;

[0110] Step S13: selecting a filter according to the high-dimensional space distance data to obtain a filter function; dividing the value range of the filter function into multiple overlapping intervals according to the high-dimensional space distance data to obtain interval division data;

[0111] Step S14: performing local clustering on the preprocessed gene expression matrix according to the interval division data to obtain a local clustering result;

[0112] Step S15: constructing a cell topology map according to the local clustering results to obtain a cell topology map.

[0113] In an embodiment of the present invention, first, the original gene expression matrix is ​​read from the sequencing instrument output file, the matrix is ​​a gene in rows, a cell in columns, and the matrix elements represent the expression of each gene in each cell. The original gene expression matrix is ​​usually stored in a format such as CSV, TSV or HDF5. After reading, the data is loaded into a DataFrame format using Python's Pandas library. Next, the original gene expression matrix is ​​preprocessed. The preprocessing steps include: using the `scanpy.pp.filter_cells` function in the Scanpy library to filter out low-quality cells (e.g., cells with too few expressed genes or too high a proportion of mitochondrial genes); using the `scanpy.pp.filter_genes` function to filter out low-expressed genes (e.g., genes expressed in very few cells); using the `scanpy.pp.normalize_total` function to normalize cell expression to adjust the total expression of each cell to the same value; using the `scanpy.pp.log1p` function to perform logarithmic transformation to convert the expression to a logarithmic scale; and finally, using the `scanpy.pp.highly_variable_genes` function to select highly variable genes, which have large expression differences between different cells and better reflect cell heterogeneity. The preprocessed gene expression matrix still has rows as genes and columns as cells. The matrix elements represent the normalized and log-transformed gene expression, and only contain highly variable genes.

[0114] The high-dimensional spatial distance between cells in the preprocessed gene expression matrix is ​​calculated using Python's sklearn library. First, the preprocessed gene expression matrix is ​​converted into NumPy array format. Then, the Euclidean distance between cells is calculated using the `sklearn.metrics.pairwise_distances` function. The input of this function is the gene expression matrix, and the output is a cell-to-cell distance matrix, where each element (i, j) represents the Euclidean distance between cell i and cell j.

[0115] According to the high-dimensional space distance data, the Gaussian kernel function is selected as the filter function, and its formula is: $K(x, y) = exp(-||xy||^2 / 2σ^2)$, where x and y represent the coordinates of two cells in the high-dimensional space (i.e., gene expression vectors), $||xy||$ represents the Euclidean distance between x and y, and σ is the bandwidth parameter, whose value is set to the 75% quantile of the inter-cell distance distribution. The value range of the Gaussian kernel function [0, 1] is divided into 10 overlapping intervals: [0, 0.2], [0.1, 0.3], [0.2, 0.4], ..., [0.8, 1.0]. The distance between each cell pair is substituted into the Gaussian kernel function to calculate the function value, and the cell pairs are assigned to the corresponding overlapping intervals according to the function value.

[0116] The cells that fall within each window in the preprocessed gene expression matrix are clustered according to the interval partitioning data. Clustering is performed using the Louvain algorithm, which is a modularity-based graph clustering algorithm. First, an undirected graph is constructed based on the cell pairs that fall within each overlapping interval, where nodes represent cells and edges represent connections between cell pairs. The weight of the edge is determined by the Gaussian kernel function value. Then, the Louvain algorithm is executed using the `leidenalg.find_partition` function in Python's leidenalg library to obtain the clustering results within each overlapping interval.

[0117] Integrate the local clustering results of all overlapping intervals to construct a cell topology map. Initialize an empty graph with nodes representing all cells. Traverse the clustering results of each overlapping interval. If two cells are in the same cluster, add an edge connecting the two cells to the graph. The weight of the edge is set to the average value of the Gaussian kernel function in all overlapping intervals containing the two cells. The final cell topology map is an undirected weighted graph that reflects the connection strength between cells.

[0118] Preferably, step S2 comprises the following steps:

[0119] Step S21: Calculate the cell division rate according to the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculate the RNA rate of each cell using the pre-processed gene expression matrix to obtain RNA rate data;

[0120] Step S22: adding node attributes to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map;

[0121] Step S23: Initializing the node weights according to the enhanced cell topology map to obtain initial node weights;

[0122] Step S24: adjusting the node weights of the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain the division rate adjustment weights; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain the RNA rate adjustment edge weights;

[0123] Step S25: integrating the node weights and edge weights of the division rate adjustment weights and the RNA rate adjustment edge weights according to the enhanced cell topology map to obtain an integrated weight;

[0124] Step S26: applying an optimal path algorithm according to the enhanced cell topology map and the integrated weights to obtain optimal path data;

[0125] Step S27: constructing a dynamic trajectory map for the enhanced cell topology map according to the optimal path data to obtain a dynamic trajectory map.

[0126] As an example of the present invention, refer to Figure 2 As shown, in this example, step S2 includes:

[0127] Step S21: Calculate the cell division rate according to the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculate the RNA rate of each cell using the pre-processed gene expression matrix to obtain RNA rate data;

[0128] In an embodiment of the present invention, first, known cell cycle-related genes are extracted from the pre-processed gene expression matrix (for example, obtained from the CycleMarkers database). The average expression of each cell on these genes is calculated as a cell cycle score. The cell cycle score is calculated using the `omicsbean.cell_cycle_score` function of the OmicsBean library. Each cell is divided into four cycle stages, G1 / S / G2 / M, according to the cell cycle score. The cell division rate is calculated using cell cycle stage information and cell topology. For each cell, the proportion of its neighbor cells in the G2 / M stage is calculated, and the proportion is used as the division rate of the cell. Subsequently, the RNA rate is inferred using the pre- and post-splicing mRNA information in the pre-processed gene expression matrix. The RNA rate of each cell is calculated using the velocity library. Velocyto uses the RNA velocity model to infer the expression of spliced ​​and unspliced ​​mRNA of each gene, and calculates the RNA rate. The calculated RNA rate matrix is ​​smoothed, and the local regression smoothing method is used to reduce noise, and the smoothing window size is set to 5.

[0129] Step S22: adding node attributes to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map;

[0130] In an embodiment of the present invention, the cell division rate and RNA rate calculated in step S21 are added to the node attributes of the cell topology map. The graph data is operated using the Python NetworkX library. For each cell (i.e., each node in the graph), two new attributes are added: "division rate" and "RNA rate". The value of the "division rate" attribute is the division rate of the cell calculated in step S21, and the value of the "RNA rate" attribute is a vector representing the RNA rate of each gene in the cell. The cell topology map after adding the attributes is the enhanced cell topology map.

[0131] Step S23: Initializing the node weights according to the enhanced cell topology map to obtain initial node weights;

[0132] In the embodiment of the present invention, the node weights are initialized according to the enhanced cell topology map, and the initial weights of all nodes are set to 1. This means that in the initial state, all cells have the same weight in the dynamic trajectory inference.

[0133] Step S24: adjusting the node weights of the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain the division rate adjustment weights; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain the RNA rate adjustment edge weights;

[0134] In an embodiment of the present invention, first, the average division rate of the neighboring nodes of each node in the enhanced cell topology is calculated. Then, the minimum value among the average division rates of all nodes is set to the minimum division rate. Next, the initial node weight is adjusted. For each node, its division rate is divided by the minimum division rate, and the result is used as the division rate adjustment weight of the node. For edge weight adjustment, the average RNA rate vector of each node is first calculated. Then, for each edge in the graph, the Euclidean distance between the average RNA rate vectors of the two nodes connecting the edge is calculated, and the distance is used as the RNA rate difference of the edge. The minimum value among the RNA rate differences of all edges is set to the minimum RNA rate difference. Finally, the directional weight adjustment factor of each edge is calculated. For each edge, the RNA rate difference of the two nodes connecting the edge is calculated and divided by the minimum RNA rate difference, and the result is used as the directional weight adjustment factor of the edge. The directional weight adjustment factor of the edge is multiplied by the RNA rate difference of the edge to obtain the RNA rate adjustment edge weight.

[0135] Step S25: integrating the node weights and edge weights of the division rate adjustment weights and the RNA rate adjustment edge weights according to the enhanced cell topology map to obtain an integrated weight;

[0136] In the embodiment of the present invention, the split rate adjustment weight and the RNA rate adjustment edge weight are integrated into an integrated weight. For each edge, the product of the split rate adjustment weights of the two nodes connecting the edge is multiplied by the RNA rate adjustment edge weight of the edge to obtain the integrated weight of the edge.

[0137] Step S26: applying an optimal path algorithm according to the enhanced cell topology map and the integrated weights to obtain optimal path data;

[0138] In an embodiment of the present invention, the Dijkstra algorithm is used to apply the optimal path algorithm on the enhanced cell topology map to calculate the optimal path between cells. The Dijkstra algorithm is an algorithm for calculating the single-source shortest path in a graph. The integrated weight is used as the weight of the edge to calculate the shortest path between all cell pairs in the graph, and the result is stored as the optimal path data.

[0139] Step S27: constructing a dynamic trajectory map for the enhanced cell topology map according to the optimal path data to obtain a dynamic trajectory map;

[0140] In the embodiment of the present invention, a dynamic trajectory graph is constructed based on the optimal path data. The dynamic trajectory graph is a directed weighted graph. For each path in the optimal path data, the edge on the path is added to the dynamic trajectory graph. The direction of the edge is consistent with the direction of the path, and the weight of the edge is the number of times the edge appears in the optimal path data.

[0141] Preferably, step S21 includes the following steps:

[0142] Step S211: extracting the cell cycle-related gene expression matrix from the preprocessed gene expression matrix to obtain a cycle gene expression matrix;

[0143] Step S212: performing single cell cycle scoring according to the cycle gene expression matrix to obtain a single cell cycle score;

[0144] Step S213: dividing the cell cycle stages according to the single cell cycle scores to obtain the cell cycle stages;

[0145] Step S214: Calculate the cell division rate according to the cell topology and the cell cycle stage to obtain the cell division rate;

[0146] Step S215: constructing a pre- and post-splicing mRNA matrix for the pre-processed gene expression matrix to obtain a pre- and post-splicing mRNA matrix;

[0147] Step S216: Infer the RNA rate based on the mRNA matrix before and after splicing to obtain the RNA rate matrix; perform rate smoothing and noise reduction on the RNA rate matrix to obtain RNA rate data.

[0148] In the embodiment of the present invention, the expression data of cell cycle related genes are extracted from the pre-treatment gene expression matrix. Using the CycleMarkers database, the gene list related to the human cell cycle is obtained, including marker genes in three stages of G1 / S phase, S phase and G2 / M phase. The row corresponding to these marker genes in the pre-treatment gene expression matrix is ​​extracted to form a cycle gene expression matrix. Each row of the cycle gene expression matrix represents a cell cycle related gene, and each column represents a cell, and the matrix element represents the expression amount of the gene in the cell.

[0149] The cell cycle score of each cell is calculated using the `omicsbean.cell_cycle_score` function in the OmicsBean library. This function inputs the cycle gene expression matrix and the gene list provided by the CycleMarkers database (distinguishing the three phases of G1 / S, S, and G2 / M), and outputs the scores of the three phases of G1 / S, S, and G2 / M for each cell. These three scores are combined into a vector as the cell cycle score of the cell.

[0150] According to the single cell cycle score, each cell is divided into the corresponding cell cycle stage. Compare the scores of the three stages G1 / S, S and G2 / M for each cell, and take the stage with the highest score as the cell cycle stage to which the cell belongs. For example, if a cell has a G1 / S score of 0.2, an S score of 0.5, and a G2 / M score of 0.8, the cell is divided into the G2 / M stage.

[0151] Calculate the cell division rate using the cell topology and cell cycle stage information. For each cell, find all its neighbor cells in the cell topology. Calculate the number of cells in the G2 / M phase among the neighbor cells and divide it by the total number of neighbor cells to get the division rate of the cell. For example, if a cell has 10 neighbor cells, 3 of which are in the G2 / M phase, then the division rate of the cell is 3 / 10=0.3.

[0152] Get the pre- and post-splicing mRNA matrix from the aligned RNA-seq data. Use the velocyto tool to process the aligned BAM file to obtain the counts of pre-splicing (immature) mRNA molecules and post-splicing (mature) mRNA molecules for each gene in each cell. Store these counts in a matrix where each row represents a gene, each column represents a cell, and each element of the matrix is ​​a vector containing two values, representing the counts of pre-splicing and post-splicing mRNAs, respectively.

[0153] Use the velocyto tool to infer RNA rates based on the pre- and post-splicing mRNA matrix. Based on a kinetic model, velocyto infers the RNA rate of each gene in each cell by fitting the counts of pre- and post-splicing mRNA. Each row of the resulting RNA rate matrix represents a gene, each column represents a cell, and the matrix elements represent the RNA rate of the gene in the cell. Smooth and denoise the RNA rate matrix. Use the local regression smoothing method to smooth the RNA rate of each gene in all cells to reduce the impact of noise. The smoothing window size is set to 5, that is, the RNA rates of 5 cells around each gene are used for smoothing. The smoothed RNA rate matrix is ​​the RNA rate data.

[0154] Preferably, step S24 includes the following steps:

[0155] Step S241: extracting the average node division rate of the enhanced cell topology map to obtain the average node division rate; setting the minimum division rate according to the average node division rate to obtain the minimum division rate;

[0156] Step S242: adjusting the initial node weights according to the average node splitting rate and the minimum node splitting rate to obtain a splitting rate adjustment weight;

[0157] Step S243: extracting the node average RNA rate vector of the enhanced cell topology map to obtain the node average RNA rate vector; performing edge RNA rate difference calculation on the enhanced cell topology map according to the node average RNA rate vector to obtain edge RNA rate difference data;

[0158] Step S244: setting the minimum RNA rate difference according to the edge RNA rate difference data to obtain the minimum RNA rate difference data;

[0159] Step S245: calculating the directional weight adjustment factor of the edge of the enhanced cell topology map according to the edge RNA rate difference data to obtain the directional weight adjustment factor;

[0160] Step S246: adjusting the edge weight of the edge RNA rate difference data according to the minimum RNA rate difference data and the directional weight adjustment factor to obtain the RNA rate adjusted edge weight.

[0161] In the embodiment of the present invention, each node in the enhanced cell topology map is traversed. For each node, the average value of the splitting rates of its neighboring nodes is calculated, and the average value is used as the average splitting rate of the node. The average splitting rates of all nodes are stored in a list. The minimum value in the list is obtained and set as the minimum splitting rate. If the minimum splitting rate is 0, it is set to the minimum value of the non-zero value in the list to avoid division by zero in subsequent calculations.

[0162] Traverse each node in the enhanced cell topology. For each node, divide the average division rate of the node by the minimum division rate. The result is used as the division rate adjustment weight of the node. The division rate adjustment weights of all nodes are stored in a list, and the index of the list corresponds to the node index.

[0163] Traverse each node in the enhanced cell topology graph. For each node, calculate the average of the RNA rate vectors of its neighboring nodes, and use the average vector as the average RNA rate vector of the node. Traverse each edge in the enhanced cell topology graph. For each edge, get the average RNA rate vector of the two nodes connecting the edge. Calculate the Euclidean distance between the two vectors and use the distance as the RNA rate difference of the edge. Store the RNA rate differences of all edges in a list, and the index of the list corresponds to the edge index.

[0164] Get the RNA rate difference data for all edges and store them in a list. Calculate the minimum value in that list and set that minimum value as the minimum RNA rate difference. If the minimum value is 0, set it to the minimum of the non-zero values ​​in the list.

[0165] Traverse each edge in the enhanced cell topology graph. For each edge, divide the RNA rate difference of the edge by the minimum RNA rate difference. Use the calculated result as the directional weight adjustment factor of the edge. Store the directional weight adjustment factors of all edges in a list, and the index of the list corresponds to the edge index.

[0166] Traverse each edge in the enhanced cell topology graph. For each edge, multiply the RNA rate difference of the edge by the directional weight adjustment factor of the edge. The result is used as the RNA rate adjustment edge weight of the edge. Store the RNA rate adjustment edge weights of all edges in a list, and the index of the list corresponds to the edge index.

[0167] Preferably, step S3 comprises the following steps:

[0168] Step S31: performing feature matrix conversion on the dynamic trajectory graph to obtain a trajectory feature matrix;

[0169] Step S32: determining the number of clusters for the trajectory feature matrix to obtain cluster number data;

[0170] Step S33: Perform fuzzy C-means clustering according to the trajectory feature matrix and the cluster quantity data to obtain a fuzzy clustering result;

[0171] Step S34: construct a membership matrix based on the fuzzy clustering results to obtain a cell membership matrix.

[0172] In an embodiment of the present invention, the dynamic trajectory graph is converted into a trajectory feature matrix. First, the dynamic trajectory graph is represented as an adjacency matrix, in which the matrix element (i, j) represents the weight of the edge from cell i to cell j. If there is no edge from i to j, the value of the matrix element (i, j) is 0. Then, the adjacency matrix is ​​eigendecomposed to obtain eigenvalues ​​and eigenvectors. The eigenvectors corresponding to the first K largest eigenvalues ​​are selected to form a trajectory feature matrix. The K value is determined based on the variance explanation rate of the eigenvalue, and the minimum K value that can explain 90% of the variance is selected. Each row of the trajectory feature matrix represents a cell, each column represents a feature, and the matrix element represents the value of the cell on the feature.

[0173] The silhouette coefficient method is used to determine the optimal number of clusters. For a series of preset cluster number K values ​​(for example, from 2 to 10), the trajectory feature matrix is ​​clustered using the fuzzy C-means clustering algorithm. For each K value, the silhouette coefficient of the clustering result is calculated. The range of the silhouette coefficient is [-1, 1], and the larger the value, the better the clustering effect. The K value with the highest silhouette coefficient is selected as the optimal number of clusters.

[0174] Fuzzy C-means clustering is performed on the trajectory feature matrix using the FuzzyCMeans function in the sklearn.fuzzy module of Python. The input parameters include the trajectory feature matrix, the number of clusters K, the fuzzy coefficient m (set to 2), the maximum number of iterations (set to 1000), and the termination condition (set to 1e-5). The function outputs the fuzzy clustering results, including the membership of each cell to each cluster.

[0175] The cell membership matrix is ​​constructed based on the fuzzy clustering results. Each row of the cell membership matrix represents a cell, each column represents a cluster, and the matrix element (i, j) represents the membership of cell i to cluster j. The value range of membership is [0, 1]. The larger the value, the greater the possibility that cell i belongs to cluster j.

[0176] Preferably, step S4 comprises the following steps:

[0177] Step S41: clustering the preprocessed gene expression matrix according to the cell membership matrix to obtain a gene expression grouping list;

[0178] Step S42: performing clustering probability distribution calculation according to the gene expression grouping list, and performing gene sorting to obtain gene entropy sorting data;

[0179] Step S43: selecting a co-expression measurement method according to the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters;

[0180] Step S44: Calculate the inter-gene co-expression values ​​of the pre-processed gene expression matrix using the co-expression measurement method and parameters to obtain a gene co-expression matrix;

[0181] Step S45: constructing a co-expression network according to the gene co-expression matrix to obtain a gene co-expression network; identifying network modules of the gene co-expression network to obtain network modules;

[0182] Step S46: Calculate the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain the module information entropy ranking data;

[0183] Step S47: Screen characteristic genes according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list.

[0184] As an example of the present invention, refer to Figure 3 As shown, in this example, step S4 includes:

[0185] Step S41: clustering the preprocessed gene expression matrix according to the cell membership matrix to obtain a gene expression grouping list;

[0186] In the embodiment of the present invention, according to the cell membership matrix and the pre-processing gene expression matrix, the gene expression data are grouped according to clustering. For each cluster, the gene expression data of the cells belonging to the cluster are extracted. Specifically, for cluster k, cells with a membership greater than 0.5 in the cell membership matrix are found, and then the columns corresponding to these cells are extracted from the pre-processing gene expression matrix to form a new gene expression matrix. The gene expression matrices of all clusters are stored in a list, which is called the gene expression grouping list.

[0187] Step S42: performing clustering probability distribution calculation according to the gene expression grouping list, and performing gene sorting to obtain gene entropy sorting data;

[0188] In an embodiment of the present invention, for each gene expression matrix in the gene expression grouping list, the probability distribution of each gene in the cluster is calculated. The probability distribution is calculated using a Gaussian kernel density estimation method. First, a Gaussian kernel function is determined. Then, the bandwidth parameter is determined using the Silverman's rule of thumb method. Using the Gaussian kernel function and the bandwidth parameter, a kernel density estimation is performed on the expression value of each gene in the cluster to obtain the probability density function of the gene. The probability density function is normalized to obtain the probability distribution of the gene. Next, the information entropy of each gene is calculated. Based on the probability distribution of the gene, the information entropy of each gene is calculated. Information entropy is used to measure the degree of discreteness of gene expression distribution. Then, the gene information entropy is weighted according to the cell membership matrix. For each gene, its information entropy is multiplied by the sum of the membership of the corresponding cell in the cell membership matrix to obtain the weighted gene information entropy. Finally, the genes are sorted according to the weighted gene information entropy to obtain gene entropy sorting data, and the sorting order is information entropy from large to small.

[0189] Step S43: selecting a co-expression measurement method according to the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters;

[0190] In the present embodiment, the Pearson correlation coefficient is selected as the co-expression measurement method. The Pearson correlation coefficient is used to measure the linear correlation between two gene expression values. The parameters are set as follows: the absolute value of the correlation coefficient is greater than 0.8 and the p value is less than 0.05.

[0191] Step S44: Calculate the inter-gene co-expression values ​​of the pre-processed gene expression matrix using the co-expression measurement method and parameters to obtain a gene co-expression matrix;

[0192] In the embodiment of the present invention, the co-expression value between genes in the pre-processed gene expression matrix is ​​calculated using the Pearson correlation coefficient method. For each pair of genes, the Pearson correlation coefficient and the p value of their expression values ​​in all cells are calculated. The correlation coefficient and the p value are stored in the gene co-expression matrix. The gene co-expression matrix is ​​a symmetrical matrix, and the matrix element (i, j) represents the Pearson correlation coefficient between gene i and gene j.

[0193] Step S45: constructing a co-expression network according to the gene co-expression matrix to obtain a gene co-expression network; identifying network modules of the gene co-expression network to obtain network modules;

[0194] In an embodiment of the present invention, a gene coexpression network is constructed according to a gene coexpression matrix. A gene coexpression network is an undirected weighted graph in which nodes represent genes and edges represent the coexpression relationship between genes. For a gene pair (i, j) that satisfies the parameter conditions set in step S43 in the gene coexpression matrix, an edge connecting gene i and gene j is added to the gene coexpression network, and the weight of the edge is the absolute value of the Pearson correlation coefficient between gene i and gene j. The gene coexpression network is identified as a network module using the Louvain algorithm. The Louvain algorithm is a graph clustering algorithm based on modularity. The identified network modules are stored in a list, and each element of the list represents a module, and the module contains genes belonging to the module.

[0195] Step S46: Calculate the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain the module information entropy ranking data;

[0196] In the embodiment of the present invention, for each network module, the module information entropy is calculated. First, according to the cell membership matrix and the pre-processed gene expression matrix, the average expression of the genes in the module in each cluster is calculated. Then, the probability distribution of the module in each cluster is calculated based on the average expression. The module information entropy is calculated using the same method as in step S42. The information entropy of all modules is stored in a list, which is called the module information entropy sorting data, and the sorting order is the information entropy from large to small.

[0197] Step S47: Screening characteristic genes according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list;

[0198] In the embodiment of the present invention, characteristic gene screening is performed based on gene entropy sorting data and module information entropy sorting data. First, the genes ranked in the top 10% in the gene entropy sorting data are selected as candidate characteristic genes. Then, among the candidate characteristic genes, the genes belonging to the modules ranked in the top 10% in the module information entropy sorting data are selected as the final characteristic genes. The final characteristic genes are stored in a list, which is called a characteristic gene list.

[0199] Preferably, step S42 includes the following steps:

[0200] Step S421: determining a kernel function according to the gene expression grouping list to obtain a kernel function; determining a bandwidth according to the gene expression grouping list and the kernel function to obtain a bandwidth parameter;

[0201] Step S422: performing kernel density estimation on the gene expression grouping list according to the kernel function and the bandwidth parameter to obtain the gene clustering probability density; performing probability density normalization on the gene clustering probability density to obtain the gene clustering probability distribution data;

[0202] Step S423: Calculate gene information entropy according to gene clustering probability distribution data to obtain gene information entropy;

[0203] Step S424: adjusting the information entropy weight of the gene information entropy according to the cell membership matrix to obtain weighted gene information entropy;

[0204] Step S425: performing gene entropy sorting according to the weighted gene information entropy to obtain gene entropy sorting data.

[0205] In an embodiment of the present invention, a Gaussian kernel function is selected as the kernel function for kernel density estimation. The formula of the Gaussian kernel function is: $K(x)=(1 / (σ*sqrt(2π)))*exp(-(x^2) / (2*σ^2))$, where x is a variable and σ is a bandwidth parameter. For each gene expression matrix in the gene expression grouping list, the bandwidth parameter is determined using the Silverman's rule of thumb method. The formula for calculating the bandwidth using the Silverman's rule of thumb method is: $σ=0.9*min(std, IQR / 1.34)*n^(-1 / 5)$, where std is the standard deviation of the gene expression value, IQR is the interquartile range of the gene expression value, and n is the number of cells.

[0206] For each gene expression matrix in the gene expression grouping list, the Gaussian kernel function and the bandwidth parameter determined in step S421 are used to carry out kernel density estimation. For each gene, the expression value of this gene in all cells is used as input, and the kernel density estimation method is applied to calculate the probability density of this gene at each expression value. The probability density at all expression values ​​is stored in a vector, which is the gene clustering probability density. Then, the gene clustering probability density is normalized. The gene clustering probability density vector is divided by the sum of all its elements to obtain the gene clustering probability distribution data.

[0207] Gene information entropy is calculated based on gene clustering probability distribution data. For each gene, its information entropy is calculated using the following formula: $H(X)=-Σ(p(x)*log2(p(x)))$, where p(x) is the probability of the gene at expression value x, and log2 is the logarithm with base 2. The probability distribution of each gene is traversed, each probability value and its corresponding logarithm are multiplied and summed, and finally the negative value is taken to obtain the gene information entropy.

[0208] The information entropy weight of the gene information entropy is adjusted according to the cell membership matrix. For each gene, its information entropy is multiplied by the sum of the membership of the corresponding cell in the cell membership matrix. Specifically, for gene i, $∑(j=1toC)M(c,j)$ is calculated, where C is the number of clusters, M is the cell membership matrix, c is the cell index, and j is the cluster index. The calculated result is multiplied by the information entropy of gene i to obtain the weighted gene information entropy.

[0209] Gene entropy sorting is performed according to weighted gene information entropy. All genes are sorted from large to small according to weighted gene information entropy to obtain gene entropy sorting data. Gene entropy sorting data contains genes and their corresponding weighted gene information entropy, and are arranged in descending order according to entropy values.

[0210] Preferably, step S5 comprises the following steps:

[0211] Step S51: preparing training data according to the characteristic gene list and the preset cell marker gene database to obtain a training data set;

[0212] Step S52: training a multi-label classification model using the training data set to obtain a multi-label classification model;

[0213] Step S53: performing model evaluation and optimization on the multi-label classification model to obtain an optimized multi-label classification model;

[0214] Step S54: According to the preprocessed gene expression matrix and the characteristic gene list, the optimized multi-label classification model is used to predict the cell state to obtain the cell state prediction result;

[0215] Step S55: Generate a cell state annotation result according to the cell state prediction result to obtain a cell state annotation result.

[0216] In an embodiment of the present invention, the CellMarker database is used as a preset cell marker gene database. Marker gene data is downloaded from the CellMarker database, which contains a list of marker genes for different cell types. A training data set is prepared based on the characteristic gene list and the CellMarker database. For each cell type in the CellMarker database, a training sample is created. The feature vector of the sample consists of the expression values ​​of the genes in the characteristic gene list in the cell type. If a certain characteristic gene does not have expression information of the cell type in the CellMarker database, the expression value of the gene is set to 0. The label of the sample is a binary vector, and the length of the vector is equal to the number of cell types in the CellMarker database. If the sample belongs to a certain cell type, the corresponding label is set to 1, otherwise it is set to 0.

[0217] Train a multi-label classification model using the training dataset. Select RandomForestClassifier as the multi-label classification model. Create a model object using the RandomForestClassifier class from the scikit-learn library. Set the model parameters: n_estimators = 100, max_depth = 10, random_state = 0. Fit the model using the training dataset.

[0218] Use 5-fold cross validation to evaluate the performance of the multi-label classification model. Use the cross_val_score function in the scikit-learn library for cross validation. Use Hamming loss as the evaluation indicator. The value range of Hamming loss is [0, 1]. The smaller the value, the better the model performance. Use the grid search method to optimize the model parameters. The search parameter space is: n_estimators = [50, 100, 150], max_depth = [5, 10, 15]. Select the parameter combination with the smallest Hamming loss as the best parameter. Use the best parameters to retrain the multi-label classification model to obtain the optimized multi-label classification model.

[0219] Based on the preprocessed gene expression matrix and the list of characteristic genes, the optimized multi-label classification model is used to predict the cell state. First, the rows corresponding to the characteristic genes are extracted from the preprocessed gene expression matrix to form a new gene expression matrix. Then, the new gene expression matrix is ​​used as input and the optimized multi-label classification model is used for prediction. The model outputs the predicted label of each cell, which is a binary vector whose length is equal to the number of cell types in the CellMarker database.

[0220] Generate cell state annotation results based on cell state prediction results. For each cell, find the cell type corresponding to the label with a value greater than 0.5 in its predicted label. Use these cell types as the annotation results of the cell. If a cell does not have a label with a value greater than 0.5 in its predicted label, annotate the cell as "Unknown". Store the annotation results of all cells in a list, called cell state annotation results.

[0221] Preferably, the present invention further provides a single-cell sequencing data classification system, which is used to perform the single-cell sequencing data classification method as described above, and the single-cell sequencing data classification system comprises:

[0222] A manifold embedding module is used to obtain an original gene expression matrix; perform gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; and construct a cell topology map according to the preprocessed gene expression matrix to obtain a cell topology map;

[0223] The dynamic trajectory inference module is used to calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights according to the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference according to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights to obtain a dynamic trajectory map;

[0224] The fuzzy clustering module is used to perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix;

[0225] The characteristic gene screening module is used to calculate the clustering probability distribution according to the cell membership matrix and the preprocessed gene expression matrix, and perform gene sorting to obtain gene entropy sorting data; calculate the module information entropy according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; perform characteristic gene screening according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list;

[0226] The cell state annotation module is used to annotate the cell state according to the characteristic gene list and obtain the cell state annotation result to realize the data classification task of single-cell sequencing.

[0227] Therefore, the embodiments should be regarded as illustrative and non-restrictive from all points of view, and the scope of the present invention is limited by the appended claims rather than the above description, and it is intended that all changes falling within the meaning and scope of equivalent elements of the application documents are included in the present invention.

[0228] The above description is only a specific embodiment of the present invention, so that those skilled in the art can understand or implement the present invention. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but should conform to the widest scope consistent with the principles and novel features invented herein.

Claims

1. A data classification method for single-cell sequencing, characterized in that: The following steps are involved: Step S1: Obtain the original gene expression matrix; Performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; The cell topology map is constructed according to the preprocessed gene expression matrix to obtain the cell topology map; Step S2: Calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain cell division rate data and RNA rate data; add node rates to the cell topology map according to the cell division rate data and RNA rate data to obtain an enhanced cell topology map; According to the enhanced cell topology, node weights are adjusted and edge weights are adjusted to obtain division rate adjustment weights and RNA rate adjustment edge weights; According to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights, dynamic trajectory inference is performed to obtain a dynamic trajectory map; Step S3: Perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix; Step S4: Calculate clustering probability distribution based on the cell membership matrix and the preprocessed gene expression matrix, and sort genes to obtain gene entropy sorting data; Calculate the module information entropy based on the cell membership matrix and the preprocessed gene expression matrix to obtain the module information entropy sorting data; screen the characteristic genes based on the gene entropy sorting data and the module information entropy sorting data to obtain the characteristic gene list; Step S5: Annotate the cell state according to the characteristic gene list to obtain the cell state annotation result to realize the data classification task of single-cell sequencing.

2. The data classification method for single-cell sequencing according to claim 1, characterized in that: Step S1 includes the following steps: Step S11: obtaining an original gene expression matrix; performing gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; Step S12: performing high-dimensional space distance calculation on the distances between cells in the preprocessed gene expression matrix to obtain high-dimensional space distance data; Step S13: selecting a filter according to the high-dimensional space distance data to obtain a filter function; dividing the value range of the filter function into multiple overlapping intervals according to the high-dimensional space distance data to obtain interval division data; Step S14: performing local clustering on the preprocessed gene expression matrix according to the interval division data to obtain a local clustering result; Step S15: constructing a cell topology map according to the local clustering results to obtain a cell topology map.

3. The data classification method for single-cell sequencing according to claim 1, characterized in that: Step S2 includes the following steps: Step S21: Calculate the cell division rate according to the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculate the RNA rate of each cell using the pre-processed gene expression matrix to obtain RNA rate data; Step S22: adding node attributes to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; Step S23: Initializing the node weights according to the enhanced cell topology map to obtain initial node weights; Step S24: adjusting the node weights of the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain the division rate adjustment weights; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain the RNA rate adjustment edge weights; Step S25: integrating the node weights and edge weights of the division rate adjustment weights and the RNA rate adjustment edge weights according to the enhanced cell topology map to obtain an integrated weight; Step S26: applying an optimal path algorithm according to the enhanced cell topology map and the integrated weights to obtain optimal path data; Step S27: constructing a dynamic trajectory map for the enhanced cell topology map according to the optimal path data to obtain a dynamic trajectory map.

4. The data classification method for single-cell sequencing according to claim 3, characterized in that: Step S21 includes the following steps: Step S211: extracting the cell cycle-related gene expression matrix from the preprocessed gene expression matrix to obtain a cycle gene expression matrix; Step S212: performing single cell cycle scoring according to the cycle gene expression matrix to obtain a single cell cycle score; Step S213: dividing the cell cycle stages according to the single cell cycle scores to obtain the cell cycle stages; Step S214: Calculate the cell division rate according to the cell topology and the cell cycle stage to obtain the cell division rate; Step S215: constructing a pre- and post-splicing mRNA matrix for the pre-processed gene expression matrix to obtain a pre- and post-splicing mRNA matrix; Step S216: Infer the RNA rate based on the mRNA matrix before and after splicing to obtain the RNA rate matrix; perform rate smoothing and noise reduction on the RNA rate matrix to obtain RNA rate data.

5. The data classification method for single-cell sequencing according to claim 3, characterized in that: Step S24 includes the following steps: Step S241: extracting the average node division rate of the enhanced cell topology map to obtain the average node division rate; setting the minimum division rate according to the average node division rate to obtain the minimum division rate; Step S242: adjusting the initial node weights according to the average node splitting rate and the minimum node splitting rate to obtain a splitting rate adjustment weight; Step S243: extracting the node average RNA rate vector of the enhanced cell topology map to obtain the node average RNA rate vector; performing edge RNA rate difference calculation on the enhanced cell topology map according to the node average RNA rate vector to obtain edge RNA rate difference data; Step S244: setting the minimum RNA rate difference according to the edge RNA rate difference data to obtain the minimum RNA rate difference data; Step S245: calculating the directional weight adjustment factor of the edge of the enhanced cell topology map according to the edge RNA rate difference data to obtain the directional weight adjustment factor; Step S246: adjusting the edge weight of the edge RNA rate difference data according to the minimum RNA rate difference data and the directional weight adjustment factor to obtain the RNA rate adjusted edge weight.

6. The data classification method for single-cell sequencing according to claim 1, characterized in that: Step S3 includes the following steps: Step S31: performing feature matrix conversion on the dynamic trajectory graph to obtain a trajectory feature matrix; Step S32: determining the number of clusters for the trajectory feature matrix to obtain cluster number data; Step S33: Perform fuzzy C-means clustering according to the trajectory feature matrix and the cluster quantity data to obtain a fuzzy clustering result; Step S34: constructing a membership matrix according to the fuzzy clustering results to obtain a cell membership matrix.

7. The data classification method for single-cell sequencing according to claim 1, characterized in that: Step S4 includes the following steps: Step S41: clustering the preprocessed gene expression matrix according to the cell membership matrix to obtain a gene expression grouping list; Step S42: performing clustering probability distribution calculation according to the gene expression grouping list, and performing gene sorting to obtain gene entropy sorting data; Step S43: selecting a co-expression measurement method according to the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters; Step S44: Calculate the inter-gene co-expression values ​​of the pre-processed gene expression matrix using the co-expression measurement method and parameters to obtain a gene co-expression matrix; Step S45: constructing a co-expression network according to the gene co-expression matrix to obtain a gene co-expression network; identifying network modules of the gene co-expression network to obtain network modules; Step S46: Calculate the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain the module information entropy ranking data; Step S47: Screen characteristic genes according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list.

8. The data classification method for single-cell sequencing according to claim 7, characterized in that: Step S42 includes the following steps: Step S421: determining a kernel function according to the gene expression grouping list to obtain a kernel function; determining a bandwidth according to the gene expression grouping list and the kernel function to obtain a bandwidth parameter; Step S422: performing kernel density estimation on the gene expression grouping list according to the kernel function and the bandwidth parameter to obtain the gene clustering probability density; performing probability density normalization on the gene clustering probability density to obtain the gene clustering probability distribution data; Step S423: Calculate gene information entropy according to gene clustering probability distribution data to obtain gene information entropy; Step S424: adjusting the information entropy weight of the gene information entropy according to the cell membership matrix to obtain weighted gene information entropy; Step S425: performing gene entropy sorting according to the weighted gene information entropy to obtain gene entropy sorting data.

9. The data classification method for single-cell sequencing according to claim 1, characterized in that: Step S5 includes the following steps: Step S51: preparing training data according to the characteristic gene list and the preset cell marker gene database to obtain a training data set; Step S52: training a multi-label classification model using the training data set to obtain a multi-label classification model; Step S53: performing model evaluation and optimization on the multi-label classification model to obtain an optimized multi-label classification model; Step S54: According to the preprocessed gene expression matrix and the characteristic gene list, the optimized multi-label classification model is used to predict the cell state to obtain the cell state prediction result; Step S55: Generate a cell state annotation result according to the cell state prediction result to obtain a cell state annotation result.

10. A data classification system for single-cell sequencing, characterized in that: Used to perform the data classification method for single-cell sequencing as claimed in claim 1, the data classification system for single-cell sequencing comprises: A manifold embedding module is used to obtain an original gene expression matrix; perform gene data preprocessing on the original gene expression matrix to obtain a preprocessed gene expression matrix; and construct a cell topology map according to the preprocessed gene expression matrix to obtain a cell topology map; The dynamic trajectory inference module is used to calculate the cell division rate and RNA rate according to the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map according to the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights according to the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference according to the enhanced cell topology map, the division rate adjustment weights and the RNA rate adjustment edge weights to obtain a dynamic trajectory map; The fuzzy clustering module is used to perform fuzzy clustering analysis on the dynamic trajectory graph and construct a membership matrix to obtain a cell membership matrix; The characteristic gene screening module is used to calculate the clustering probability distribution according to the cell membership matrix and the preprocessed gene expression matrix, and perform gene sorting to obtain gene entropy sorting data; calculate the module information entropy according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; perform characteristic gene screening according to the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list; The cell state annotation module is used to annotate the cell state according to the characteristic gene list and obtain the cell state annotation result to realize the data classification task of single-cell sequencing.

Images