A data classification method and system for single-cell sequencing
By constructing cell topology maps and performing dynamic trajectory inference, combined with fuzzy clustering and characteristic gene screening, the problem of difficult to capture dynamic changes in cell state and parameter sensitivity in existing methods is solved, and more refined cell state annotation and heterogeneity assessment are achieved.
Patent Information
- Application Number
- CN202411638407.7
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2024-11-17
- Publication Date
- 2025-08-22
- Estimated Expiration
- 2044-11-17
AI Technical Summary
Existing single-cell sequencing data analysis methods are difficult to capture the dynamic changes in cell state, the clustering results are sensitive to algorithm parameters, difficult to deal with complex cell population structure, and lack quantitative assessment of cell heterogeneity.
By constructing a cell topology map, cell division rate and RNA rate are calculated, dynamic trajectory inference is performed, and cell state annotation is achieved by combining fuzzy clustering and characteristic gene screening.
Accurately reflect the connection relationship and dynamic changes between cells, reduce sensitivity to parameters, provide quantitative evaluation of cellular heterogeneity, and improve the accuracy of cell state annotation.
Smart Images

Figure CN119724350B_ABST
Abstract
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 analysis of single-cell sequencing data primarily relied on traditional clustering methods, such as k-means and hierarchical clustering. These methods group cells based on similarities in gene expression, but they often assume that cell populations are static and ignore the dynamic changes in cell state. With advances in single-cell sequencing technology and the increase in data volume, new analysis methods have emerged, such as graph-based clustering, methods based on dimensionality reduction and visualization, and trajectory inference. However, existing methods are insufficient for analyzing cellular heterogeneity.
[0003] Manifestations of insufficient analysis of cellular heterogeneity:
[0004] Difficulty capturing 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 often sensitive to the choice of algorithm parameters, such as the k value in k-means. Different parameter choices can lead to different clustering results, affecting the understanding of cellular 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 data classification method and system for single-cell sequencing to solve at least one of the above technical problems.
[0009] To achieve the above objectives, 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 based on the preprocessed gene expression matrix to obtain a cell topology map;
[0011] Step S2: Calculate the cell division rate and RNA rate based on the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map based on the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights based on the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference based on the enhanced cell topology map, division rate adjustment weights, and 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: Calculate cluster probability distribution based on the cell membership matrix and the preprocessed gene expression matrix, and perform gene sorting to obtain gene entropy sorting data; calculate module information entropy based on the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; screen characteristic genes based on the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list;
[0014] Step S5: Perform cell state annotation based on the feature gene list to obtain cell state annotation results to implement data classification for single-cell sequencing.
[0015] The present invention effectively removes noise data, highlights highly variable genes, and constructs a topology map based on local clustering results, more accurately reflecting the connection relationship between cells and the local structure in high-dimensional space, providing a more reliable and refined basis for subsequent analysis. By calculating and integrating cell division rate and RNA rate information, an enhanced cell topology map is constructed, and the dynamic trajectory map is inferred using integrated weights and the optimal path algorithm, 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 characteristic 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 characteristic genes that contribute the most to cell state annotation, improving the accuracy and biological significance of cell state annotation. By utilizing a pre-set database of cell marker genes and screened characteristic genes, training and optimizing a multi-label classification model, and predicting and annotating cell states, the automated classification of single-cell sequencing data is achieved, providing more comprehensive and accurate cell state annotation results, facilitating 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 dynamic graphs and analyzing their changing patterns, this method can better capture the dynamic changes in cell states, reduce sensitivity to parameters, handle complex cell population structures, and provide quantitative cell heterogeneity assessment indicators, thereby overcoming the shortcomings of existing methods for cell heterogeneity analysis and deepening the understanding of cell heterogeneity.
[0016] Preferably, step S1 includes 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: performing filter selection 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 based on the local clustering results to obtain a cell topology map.
[0022] The present invention can remove noise and low-quality data, such as low-expressed 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 differences 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 most to cellular heterogeneity, making the constructed cell topology map 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 maps, allowing the topology map to accurately reflect the relationship between cells. Selecting an appropriate filter function and dividing its value range into multiple overlapping intervals can effectively capture the local neighboring relationships 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 rigid partitioning. Local clustering is performed on each overlapping interval to identify the local structure between cells at different scales. This helps to construct a more refined cell topology map that can reflect the multi-level relationship between cells, rather than just based on global similarity. Constructing a cell topology based on local clustering results can more accurately reflect the connectivity between cells and preserve the local structural information of cells in high-dimensional space. Compared to constructing a topology directly based on a distance matrix, this method can better handle the noise and complex structures in high-dimensional data, making the constructed topology more robust and more representative of the actual cellular relationships.
[0023] Preferably, step S2 includes the following steps:
[0024] Step S21: Calculating the cell division rate based on the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculating 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 node weights according to the enhanced cell topology map to obtain initial node weights;
[0027] Step S24: adjusting the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain a division rate adjustment weight; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain an RNA rate adjustment edge weight;
[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 based on 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 expression of cell cycle-related genes and the mRNA information before and after splicing, which can more accurately capture the dynamic changes of cells. The cell division rate reflects the activity of cell proliferation, while the RNA rate reflects the dynamic change trend of gene expression. These two indicators provide important biological information for 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 subsequent node weight and edge weight adjustment, so that dynamic trajectory inference can better reflect the real process of cell development. Initializing the node weight provides a starting point for subsequent weight adjustment based on cell division rate and RNA rate. The initialization weight is 1, which ensures that all cells have the same weight in the initial state and avoids the influence of the initial value on the subsequent calculation results. Adjusting the node weight based on the cell division rate allows cells with high division rates to obtain higher weights in trajectory inference, highlighting the importance of these cells in the cell development process. At the same time, edge weights are adjusted based on RNA rate, giving higher weights to edges with larger RNA rate variations, thereby better capturing the impact of dynamic changes in gene expression on cell-state transitions. The division rate-adjusted weights and RNA rate-adjusted edge weights are integrated to comprehensively consider the impact of both cell division and dynamic changes in gene expression on cell-state transitions, making the integrated weights more reflective of the driving forces of cell development. Applying an optimal path algorithm to the enhanced cell topology map identifies the most likely paths connecting different cell states, representing the trajectories of cell development. Using the integrated weights as the basis for path calculation, the optimal path reflects the combined influence of cell division and dynamic changes in gene expression. A dynamic trajectory map is constructed based on the optimal path data, transforming static cell-cell connectivity into dynamic evolutionary trajectories. This more intuitive representation of cell-state transitions provides an important foundation 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: Calculating the cell division rate based on 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: RNA rate inference is performed based on the mRNA matrix before and after splicing to obtain an RNA rate matrix; rate smoothing and noise reduction are performed on the RNA rate matrix to obtain RNA rate data.
[0039] The present invention focuses on genes closely related to the cell cycle process by extracting a cell cycle-related gene expression matrix, removing the interference of other gene expression information, making the cell cycle score more accurate and more able to 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 simply relying 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 accurately describe the proliferation state of the cell. This is crucial for understanding cell heterogeneity and developmental trajectory because cells in different cycle stages have different transcriptome characteristics and biological functions. Calculating cell division rate in combination with cell topology and cell cycle stage information takes into account the relationship between the cell and its neighboring cells, and can more accurately reflect the local proliferation activity of the cell. For example, even if a cell's cycle score is high, if its neighboring cells are all in the resting phase, its actual division rate may not be high. Constructing a pre- and post-splicing mRNA matrix provides key data for RNA rate inference. Pre-splicing mRNA represents newly synthesized transcripts, while post-splicing mRNA represents mature transcripts. By comparing the ratios of the two, we can infer dynamic trends in gene expression. Inferring RNA rates based on the pre- and post-splicing mRNA matrix can capture dynamic changes in gene expression. By smoothing and noise reduction to remove the effects of technical errors, RNA rate data is more reliable and better reflects 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 from 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 based on the node average RNA rate vector to obtain edge RNA rate difference data;
[0044] Step S244: setting a 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 edge directional weight adjustment factor of the enhanced cell topology map according to the edge RNA rate difference data to obtain the directional weight adjustment factor;
[0046] Step S246: performing edge weight adjustment on the edge RNA rate difference data according to the minimum RNA rate difference data and the directional weight adjustment factor to obtain an 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 impact of individual cell division rate outliers 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. Therefore, 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 transition. 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 calculation of the directional weight adjustment factor for edges takes into account the relative magnitude 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. Edge weights are adjusted based on the minimum RNA rate difference and the directional weight adjustment factor, giving edges with large RNA rate changes a higher weight. This emphasizes the role of dynamic changes in gene expression in driving cell state transitions, making the constructed dynamic trajectory diagram more reflective of the actual process of cellular development.
[0048] Preferably, step S3 includes 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: performing fuzzy C-means clustering according to the trajectory feature matrix and the cluster quantity data to obtain the fuzzy clustering result;
[0052] Step S34: constructing a membership matrix based on the fuzzy clustering results to obtain a cell membership matrix.
[0053] The present invention converts the dynamic trajectory graph into a trajectory feature matrix, converting the graph structure information into a numerical matrix to facilitate subsequent cluster analysis. Eigendecomposition extracts the key features of the dynamic trajectory graph, reduces the dimensionality of the data, and retains key information, making the clustering results more robust and easier to interpret. Using the silhouette coefficient method to determine the number of clusters avoids the bias caused by subjective selection of the number of clusters, making the clustering results more objective and reasonable. The silhouette coefficient can effectively assess the quality of the clustering results. Selecting the K value with the highest silhouette coefficient can find the optimal number of clusters, thereby better identifying cell subpopulations. Clustering the trajectory feature matrix using the fuzzy C-means clustering algorithm can effectively handle the transition states that may exist between different cell states. Fuzzy clustering allows a cell to belong to multiple clusters simultaneously and be assigned different membership degrees, which better reflects the continuity and dynamics of cell states than hard clustering. Constructing a cell membership matrix quantifies 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 a cell belongs can be determined based on its membership degree.
[0054] Preferably, step S4 includes the following steps:
[0055] Step S41: performing clustering on 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 based on the gene expression grouping list and performing gene sorting to obtain gene entropy sorting data;
[0057] Step S43: selecting a co-expression measurement method based on the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters;
[0058] Step S44: Calculating 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; performing network module identification on the gene co-expression network to obtain a network module;
[0060] Step S46: Calculating the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy ranking data;
[0061] 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.
[0062] By grouping the gene expression matrix into clusters, the present invention can analyze gene expression patterns for each cell subpopulation separately, thereby identifying genes specifically expressed in different subpopulations. This helps understand the differences between different cell types or states and provides a foundation for subsequent signature gene screening. Calculating the gene cluster probability distribution and gene information entropy can quantify the specificity of gene expression in different cell subpopulations. The higher the information entropy, the greater the difference in gene expression between different subpopulations, and the more likely it is to be a signature gene. Weighted gene information entropy further takes into account cell affiliation, making the selection of signature genes more accurate. Selecting appropriate co-expression measurement methods and parameters can accurately assess the co-expression relationships between genes. The Pearson correlation coefficient is a commonly used co-expression measurement method that can capture linear relationships between gene expression. Appropriate parameter settings can effectively filter out weakly or falsely correlated gene pairs. Calculating inter-gene co-expression values and constructing a gene co-expression matrix provide the data foundation for subsequent construction of a gene co-expression network. The gene co-expression matrix reflects the similarity of gene expression patterns and can be used to identify functionally related gene modules. Constructing a gene co-expression network and identifying network modules can reveal gene interactions and functional modules. Modules in co-expression networks are typically composed of functionally related genes and can reflect biological pathways or regulatory mechanisms within cells. Calculating the information entropy of network modules can assess the specificity of modules across different cell subpopulations. Higher information entropy indicates greater variability in module expression across subpopulations, making it more likely to be a functional module associated with a specific cell type or state. Combining gene entropy ranking data with module entropy ranking data for signature gene screening can more effectively identify signature genes associated with cell state annotations. This approach comprehensively considers both individual and module specificity of gene expression, resulting in more representative and biologically meaningful signature genes.
[0063] Preferably, step S42 includes the following steps:
[0064] Step S421: determining a kernel function based on the gene expression group list to obtain a kernel function; determining a bandwidth based on the gene expression group 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 bandwidth parameters to obtain gene cluster probability density; performing probability density normalization on the gene cluster probability density to obtain gene cluster probability distribution data;
[0066] Step S423: Calculate gene information entropy based on gene clustering probability distribution data to obtain gene information entropy;
[0067] Step S424: adjusting the gene information entropy weight 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 effectively performs kernel density estimation by selecting a Gaussian kernel function and using Silverman's rule of thumb to determine the bandwidth. The Gaussian kernel function has excellent mathematical properties, while the Silverman's rule of thumb method can adaptively select the bandwidth based on the distribution characteristics of the data, making the kernel density estimation results more accurate. Using the kernel density estimation method to calculate the probability density of gene clusters can more accurately describe the distribution of gene expression in different cell subpopulations. Compared to simple histogram methods, kernel density estimation can better handle noise and fluctuations in the data and provide a smoother probability density curve. Probability density normalization ensures the validity of the probability distribution. Calculating gene information entropy can quantify the degree of dispersion of gene expression in different cell subpopulations. Higher information entropy indicates greater differences in gene expression between different subpopulations, and is more likely to be a characteristic gene. Using a cell membership matrix to weight gene information entropy can more accurately reflect the importance of genes in different cell subpopulations. Considering cell membership can avoid the information loss caused by fuzzy clustering and highlight genes that are highly expressed and have high information entropy in specific subpopulations. Sorting genes based on weighted gene information entropy can prioritize genes with high information entropy and high expression in specific cell subpopulations. This helps in subsequent feature gene screening, prioritizing genes that contribute most to cell state annotation.
[0070] Preferably, step S5 includes the following steps:
[0071] Step S51: preparing training data based on 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: Based on 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 based on the cell state prediction result to obtain a cell state annotation result.
[0076] The present invention prepares a training dataset using a list of characteristic genes and a preset cell marker gene database, providing high-quality annotated data for training a multi-label classification model. Reliable cell type marker gene information can be obtained using the CellMarker database. Combined with the selected characteristic genes, an information-rich and biologically meaningful training dataset can be constructed. The training dataset is used to train a multi-label classification model, enabling the model to learn the association between characteristic gene expression patterns and cell types. The RandomForestClassifier algorithm has good generalization and robustness, making it suitable for processing high-dimensional data and multi-label classification problems. Model evaluation and optimization can improve the predictive performance of the multi-label classification model. Using cross-validation and Hamming loss to evaluate the model's generalization ability, and using grid search to find the optimal model parameters, overfitting can be effectively avoided and the model's prediction accuracy can be improved. Using the optimized multi-label classification model to predict cell states on a preprocessed gene expression matrix, the cell type or state to which each cell belongs can be predicted based on the expression pattern of the characteristic genes. This provides quantitative prediction results for subsequent cell state annotation. Cell state annotation results are generated based on the cell state prediction results, and the prediction results are converted into easy-to-understand cell type or state labels. This allows researchers to easily classify and interpret single-cell data and gain insights into cellular heterogeneity and functional diversity.
[0077] Preferably, the present invention further provides a single-cell sequencing data classification system for executing the single-cell sequencing data classification method described above, the single-cell sequencing data classification system comprising:
[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 based on the preprocessed gene expression matrix to obtain a cell topology map.
[0079] A dynamic trajectory inference module is used to calculate the cell division rate and RNA rate based on the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map based on the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights based on the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference based on the enhanced cell topology map, division rate adjustment weights, and 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 the membership matrix to obtain the cell membership matrix;
[0081] The characteristic gene screening module is used to calculate the clustering probability distribution based on 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 based on the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; and screen characteristic genes based on 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 cell states based on the feature gene list and obtain cell state annotation results to implement data classification operations for single-cell sequencing. BRIEF DESCRIPTION OF THE DRAWINGS
[0083] Figure 1 A schematic diagram of the steps of a single-cell sequencing data classification method;
[0084] Figure 2 for Figure 1 Detailed implementation steps of step S2 in FIG.
[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 and in conjunction with the embodiments. DETAILED DESCRIPTION
[0087] The following is a clear and complete description of the technical method of the present invention in conjunction with the accompanying drawings. It is obvious that the embodiments described are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative efforts are within the scope of protection of the present invention.
[0088] In addition, the accompanying drawings are merely schematic illustrations of the present invention and are not necessarily drawn to scale. Identical reference numerals in the figures denote identical or similar parts, and thus repetitive descriptions thereof will be omitted. Some of the block diagrams shown in the accompanying drawings are functional entities that do not necessarily correspond to physically or logically separate entities. These functional entities may be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor and / or microcontroller approaches.
[0089] It should be understood that although the terms "first," "second," and the like may be used herein to describe various elements, these elements should not be limited by these terms. These terms are used solely to distinguish one element from another. For example, a first element may be referred to as a second element, and similarly, a second element may be referred to as a first element, without departing from the scope of the exemplary embodiments. The term "and / or" as 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 based on the preprocessed gene expression matrix to obtain a cell topology map;
[0092] Step S2: Calculate the cell division rate and RNA rate based on the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map based on the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights based on the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference based on the enhanced cell topology map, division rate adjustment weights, and 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: Calculate cluster probability distribution based on the cell membership matrix and the preprocessed gene expression matrix, and perform gene sorting to obtain gene entropy sorting data; calculate module information entropy based on the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; screen characteristic genes based on the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic gene list;
[0095] Step S5: Perform cell state annotation based on the feature gene list to obtain cell state annotation results to implement data classification for single-cell sequencing.
[0096] In the embodiment of the present invention, reference Figure 1 FIG. 1 is a flow chart showing 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 based on 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 selecting highly variable genes. Then, based on the preprocessed gene expression matrix, the Euclidean distance between cells is calculated, the Gaussian kernel function is selected as the 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 based on the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map based on the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights based on the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference based on the enhanced cell topology map, division rate adjustment weights, and RNA rate adjustment edge weights to obtain a dynamic trajectory map;
[0100] In an embodiment of the present invention, based on the cell topology map and the pre-processed gene expression matrix, the cell cycle-related gene expression matrix is extracted, the cell cycle score is calculated, and the 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, the RNA rate is calculated using velocyto, and smoothing and noise reduction processing is performed. The cell division rate and 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 based on 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 edge weight are integrated, the optimal path is calculated using the Dijkstra algorithm, and a dynamic trajectory map is constructed based on 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 this embodiment of the present invention, the dynamic trajectory graph is converted into an adjacency matrix, eigendecomposition is performed, and the top K eigenvectors (K value is determined by the variance explained) are selected to form the trajectory feature matrix. The silhouette coefficient method is used to determine the optimal number of clusters K. The trajectory feature matrix is then clustered using the fuzzy C-means clustering algorithm, with the fuzzy coefficient m set to 2. Finally, a cell membership matrix is constructed based on the clustering results, with the matrix elements representing the cell's membership to each cluster.
[0103] Step S4: Calculate cluster probability distribution based on the cell membership matrix and the preprocessed gene expression matrix, and perform gene sorting to obtain gene entropy sorting data; calculate module information entropy based on the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; screen characteristic genes based on the gene entropy sorting data and the module information entropy sorting data to obtain a characteristic 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 group, kernel density estimation is performed using the bandwidth determined by the Gaussian kernel function and the Silverman's rule of thumb method, the gene cluster probability distribution is calculated, and the gene information entropy is calculated. 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 modules are identified using the Louvain algorithm. The module information entropy of each module is calculated, and finally, characteristic genes are screened based on the gene entropy sorting data and the module information entropy sorting data.
[0105] 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;
[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 characteristic 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 includes 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: performing filter selection 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 based on the local clustering results to obtain a cell topology map.
[0113] In the present embodiment, first, the original gene expression matrix is read from the sequencer output file. The matrix is represented by rows as genes and columns as cells, and the matrix elements represent the expression level 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 the DataFrame format using the Python 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 expressing too few genes or a high 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 the same total expression value for each cell; using the `scanpy.pp.log1p` function to perform a logarithmic transformation to convert expression to a logarithmic scale; and finally, using the `scanpy.pp.highly_variable_genes` function to select highly variable genes, whose expression levels vary significantly between cells and better reflect cellular heterogeneity. The preprocessed gene expression matrix still has rows as genes and columns as cells, with the matrix elements representing the normalized and log-transformed gene expression levels and containing only highly variable genes.
[0114] The Python sklearn library was used to calculate high-dimensional spatial distances between cells in a preprocessed gene expression matrix. First, the preprocessed gene expression matrix was converted to a NumPy array format. Then, the `sklearn.metrics.pairwise_distances` function was used to calculate the Euclidean distances between cells. This function takes the gene expression matrix as input and outputs a cell-to-cell distance matrix, where each element (i, j) represents the Euclidean distance between cell i and cell j.
[0115] Based on the high-dimensional spatial distance data, a Gaussian kernel function was selected as the filter function. Its formula is: $K(x, y) = exp(-||xy||^2 / 2σ^2)$, where x and y represent the coordinates of the two cells in 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 Gaussian kernel function's range [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 interval based on the function value.
[0116] Based on the interval partitioning data, cells falling within each window in the preprocessed gene expression matrix are clustered. Clustering is performed using the Louvain algorithm, a modularity-based graph clustering algorithm. First, an undirected graph is constructed based on the cell pairs falling within each overlapping interval, where nodes represent cells and edges represent connections between cell pairs. Edge weights are determined by the value of a Gaussian kernel function. The Louvain algorithm is then executed using the `leidenalg.find_partition` function in the Python library leidenalg to obtain clustering results within each overlapping interval.
[0117] The cell topology map is constructed by integrating the local clustering results of all overlapping intervals. Initially, an empty graph is created, with nodes representing all cells. The clustering results for each overlapping interval are iterated over. If two cells are in the same cluster, an edge is added to the graph connecting them. The edge weight is set to the average of the Gaussian kernel values across all overlapping intervals containing the two cells. The resulting cell topology map is an undirected weighted graph that reflects the strength of connections between cells.
[0118] Preferably, step S2 includes the following steps:
[0119] Step S21: Calculating the cell division rate based on the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculating 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 node weights according to the enhanced cell topology map to obtain initial node weights;
[0122] Step S24: adjusting the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain a division rate adjustment weight; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain an RNA rate adjustment edge weight;
[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 based on 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: Calculating the cell division rate based on the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculating 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 the 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 neighboring cells in the G2 / M stage is calculated, and this 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 velocyto 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 Python's 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 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 initial node weights based on the cell division rate according to the enhanced cell topology map to obtain a division rate adjustment weight; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain an RNA rate adjustment edge weight;
[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 graph is calculated. Then, the minimum value among the average division rates of all nodes is set as the minimum division rate. Next, the initial node weights are 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 as 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 this embodiment of the present invention, the split rate adjustment weight and the RNA rate adjustment edge weight are integrated into a single integrated weight. For each edge, the integrated weight of the edge is obtained by multiplying the product of the split rate adjustment weights of the two nodes connecting the edge by the RNA rate adjustment edge weight of the edge.
[0137] Step S26: applying an optimal path algorithm based on the enhanced cell topology map and the integrated weights to obtain optimal path data;
[0138] In this embodiment of the present invention, the Dijkstra algorithm is applied to the enhanced cell topology map to calculate the optimal path between cells. The Dijkstra algorithm is an algorithm for calculating the shortest path from a single source in a graph. Using the integrated weights as edge weights, the shortest paths between all pairs of cells in the graph are calculated and the results are stored as 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 this 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, edges along the path are added to the dynamic trajectory graph. The direction of the edge is consistent with the path direction, 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: Calculating the cell division rate based on 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: RNA rate inference is performed based on the mRNA matrix before and after splicing to obtain an RNA rate matrix; rate smoothing and noise reduction are performed 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, comprising 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 for each cell is calculated using the `omicsbean.cell_cycle_score` function in the OmicsBean library. This function takes as input a cycle gene expression matrix and a list of genes from the CycleMarkers database (which distinguishes between G1 / S, S, and G2 / M phases). It outputs a score for each cell in the G1 / S, S, and G2 / M phases. These three scores are combined into a vector, which serves as the cell cycle score for the cell.
[0150] Each cell is assigned to the corresponding cell cycle phase based on the single cell cycle score. The scores for each cell's G1 / S, S, and G2 / M phases are compared, and the phase with the highest score is used as the cell cycle phase 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 assigned to the G2 / M phase.
[0151] Calculate cell division rates using the cell topology and cell cycle phase information. For each cell, find all of its neighboring cells in the cell topology. Count the number of neighboring cells in the G2 / M phase and divide this by the total number of neighboring cells to calculate the cell's division rate. For example, if a cell has 10 neighbors, 3 of which are in the G2 / M phase, then the cell's division rate is 3 / 10 = 0.3.
[0152] Obtain a matrix of pre- and post-splicing mRNAs from the aligned RNA-seq data. Use the velocyto tool to process the aligned BAM files and obtain counts of pre- and post-splicing (mature) mRNA molecules for each gene in each cell. These counts are stored 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, one for each pre- and post-splicing mRNA count.
[0153] The velocyto tool was used 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. The RNA rate matrix was smoothed and denoised. The local regression smoothing method was used to smooth the RNA rate of each gene in all cells to reduce the impact of noise. The smoothing window size was set to 5, that is, the RNA rates of 5 cells around each gene were 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 from 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 based on the node average RNA rate vector to obtain edge RNA rate difference data;
[0158] Step S244: setting a 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 edge directional weight adjustment factor of the enhanced cell topology map according to the edge RNA rate difference data to obtain the directional weight adjustment factor;
[0160] Step S246: performing edge weight adjustment on the edge RNA rate difference data according to the minimum RNA rate difference data and the directional weight adjustment factor to obtain an RNA rate adjusted edge weight.
[0161] In one embodiment of the present invention, each node in the enhanced cell topology is traversed. For each node, the average of the splitting rates of its neighboring nodes is calculated and 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 of the non-zero values in the list to avoid division by zero in subsequent calculations.
[0162] Traverse each node in the augmented cell topology. For each node, divide its average division rate by its minimum division rate. Use the result as the division rate adjustment weight for that node. Store all node division rate adjustment weights in a list, with the index of the list corresponding to the node index.
[0163] Traverse each node in the augmented cell map. For each node, calculate the average of the RNA rate vectors of its neighboring nodes and use this average vector as the average RNA rate vector for that node. Traverse each edge in the augmented cell map. For each edge, obtain the average RNA rate vector of the two nodes connecting that edge. Calculate the Euclidean distance between these two vectors and use this distance as the RNA rate difference for that edge. Store the RNA rate differences for all edges in a list, with the list index corresponding 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 this list and set it 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 for that edge by the minimum RNA rate difference. Use the result as the directional weight adjustment factor for that edge. Store the directional weight adjustment factors for all edges in a list, with the index of the list corresponding to the edge index.
[0166] Traverse each edge in the enhanced cell topology graph. For each edge, multiply the RNA rate difference for that edge by the directional weight adjustment factor for that edge. Use the result as the RNA rate-adjusted edge weight for that edge. Store the RNA rate-adjusted edge weights for all edges in a list, with the index of the list corresponding to the edge index.
[0167] Preferably, step S3 includes 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: performing fuzzy C-means clustering according to the trajectory feature matrix and the cluster quantity data to obtain the 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, where 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 the 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 optimal number of clusters is determined using the silhouette coefficient method. For a range of preset clustering values, K (e.g., 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 silhouette coefficient ranges from -1 to 1, with larger values indicating better clustering results. The K value with the highest silhouette coefficient is selected as the optimal number of clusters.
[0174] Fuzzy C-means clustering was performed on the trajectory feature matrix using the FuzzyCMeans function in the Python sklearn.fuzzy module. Input parameters included 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 in each cluster.
[0175] A cell membership matrix is constructed based on the fuzzy clustering results. Each row of the cell membership matrix represents a cell, and each column represents a cluster. The matrix element (i, j) represents the degree of membership of cell i to cluster j. Membership values range from [0 to 1], with larger values indicating a greater likelihood that cell i belongs to cluster j.
[0176] Preferably, step S4 includes the following steps:
[0177] Step S41: performing clustering on 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 based on the gene expression grouping list and performing gene sorting to obtain gene entropy sorting data;
[0179] Step S43: selecting a co-expression measurement method based on the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters;
[0180] Step S44: Calculating 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; performing network module identification on the gene co-expression network to obtain a network module;
[0182] Step S46: Calculating the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy ranking data;
[0183] 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.
[0184] As an example of the present invention, refer to Figure 3 As shown, in this example, step S4 includes:
[0185] Step S41: performing clustering on 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 based on 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 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 dispersion of the 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, with the sorting order being from large to small information entropy.
[0189] Step S43: selecting a co-expression measurement method based on the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters;
[0190] In the present embodiment, the Pearson correlation coefficient was selected as a co-expression measurement method. The Pearson correlation coefficient is used to measure the linear correlation between two gene expression values. The parameters were set as follows: the absolute value of the correlation coefficient was greater than 0.8 and the p-value was less than 0.05.
[0191] Step S44: Calculating 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 coexpression 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 coexpression matrix. The gene coexpression 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; performing network module identification on the gene co-expression network to obtain a network module;
[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 relationships between genes. For each 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, where 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, where each element of the list represents a module, and the module comprises the genes belonging to the module.
[0195] Step S46: Calculating the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy ranking data;
[0196] In the embodiment of the present invention, for each network module, the module information entropy is calculated. First, based on the cell membership matrix and the preprocessed gene expression matrix, the average expression of the genes in the module in each cluster is calculated. Then, based on the average expression, the probability distribution of the module in each cluster is calculated. 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, called module information entropy sorting data, and the sorting order is from large to small information entropy.
[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 this embodiment of the present invention, signature genes are screened based on gene entropy ranking data and module information entropy ranking data. First, the top 10% of genes in the gene entropy ranking data are selected as candidate signature genes. Then, among the candidate signature genes, genes belonging to the top 10% of modules in the module information entropy ranking data are selected as final signature genes. The final signature genes are stored in a list, referred to as the signature gene list.
[0199] Preferably, step S42 includes the following steps:
[0200] Step S421: determining a kernel function based on the gene expression group list to obtain a kernel function; determining a bandwidth based on the gene expression group 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 bandwidth parameters to obtain gene cluster probability density; performing probability density normalization on the gene cluster probability density to obtain gene cluster probability distribution data;
[0202] Step S423: Calculate gene information entropy based on gene clustering probability distribution data to obtain gene information entropy;
[0203] Step S424: adjusting the gene information entropy weight 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 use step S421 are carried out kernel density estimation.For each gene, the expression value of this gene in all cells is as input, and the kernel density estimation method is used to calculate the probability density of this gene at each expression value place.The probability density at all expression value places is stored in a vector, and this vector 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 the probability distribution of gene clusters. 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 base-2 logarithm. The probability distribution of each gene is traversed, and each probability value and its corresponding logarithm are multiplied and summed. The negative value is then taken to obtain the gene information entropy.
[0208] Gene entropy weighting is performed based on the cell membership matrix. For each gene, its entropy is multiplied by the sum of the memberships of the corresponding cells 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 entropy of gene i to obtain the weighted gene entropy.
[0209] Gene entropy sorting is performed based on weighted gene information entropy. All genes are sorted from largest to smallest according to their weighted gene information entropy to obtain gene entropy sorting data. Gene entropy sorting data contains genes and their corresponding weighted gene information entropy, sorted in descending order of entropy value.
[0210] Preferably, step S5 includes the following steps:
[0211] Step S51: preparing training data based on 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: Based on 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 based on 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 for 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 the Hamming loss as the evaluation metric. The Hamming loss range is [0, 1], with smaller values indicating better model performance. Use grid search to optimize the model parameters. The search parameter space is: n_estimators = [50, 100, 150], max_depth = [5, 10, 15]. The parameter combination that minimizes the Hamming loss is selected as the optimal parameter combination. Retrain the multi-label classification model using the optimal parameters to obtain the optimized multi-label classification model.
[0219] Based on the preprocessed gene expression matrix and the list of characteristic genes, an optimized multi-label classification model is used to predict cell states. First, rows corresponding to characteristic genes are extracted from the preprocessed gene expression matrix to construct a new gene expression matrix. Then, the optimized multi-label classification model is used to perform predictions using the new gene expression matrix as input. The model outputs the predicted label for each cell, a binary vector whose length equals the number of cell types in the CellMarker database.
[0220] Generate cell state annotations based on the cell state predictions. For each cell, find the cell types corresponding to the predicted labels with a value greater than 0.5. Use these cell types as the annotations for that cell. If a cell's predicted labels do not have a value greater than 0.5, annotate that cell as "Unknown." Store all cell annotations in a list called the cell state annotations.
[0221] Preferably, the present invention further provides a single-cell sequencing data classification system for executing the single-cell sequencing data classification method described above, the single-cell sequencing data classification system comprising:
[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 based on the preprocessed gene expression matrix to obtain a cell topology map.
[0223] A dynamic trajectory inference module is used to calculate the cell division rate and RNA rate based on the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map based on the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights based on the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference based on the enhanced cell topology map, division rate adjustment weights, and 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 the membership matrix to obtain the cell membership matrix;
[0225] The characteristic gene screening module is used to calculate the clustering probability distribution based on 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 based on the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; and screen characteristic genes based on 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 cell states based on the feature gene list and obtain cell state annotation results to implement data classification operations for single-cell sequencing.
[0227] Therefore, the embodiments should be regarded in all respects as illustrative and non-restrictive, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all changes that fall within the meaning and scope of equivalent elements of the application documents are included in the present invention.
[0228] The foregoing description is intended only to provide specific embodiments of the present invention, which will enable those skilled in the art to understand and implement the present invention. Various modifications to these embodiments will be readily 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 is not intended to be limited to the embodiments shown herein, but is to be construed in the widest possible manner consistent with the principles and novel features disclosed 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 based on the preprocessed gene expression matrix to obtain the cell topology map; Step S2: Calculating the cell division rate and RNA rate based on the cell topology map and the preprocessed gene expression matrix to obtain cell division rate data and RNA rate data; adding node rates to the cell topology map based on 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; Dynamic trajectory inference is performed based on the enhanced cell topology, division rate adjustment weights, and RNA rate adjustment edge weights to obtain a dynamic trajectory graph; 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 perform gene ranking to obtain gene entropy ranking data; Calculate the module information entropy based on the cell membership matrix and the preprocessed gene expression matrix to obtain the module information entropy ranking data; screen the characteristic genes based on the gene entropy ranking data and the module information entropy ranking data to obtain the characteristic gene list; Step S5: Perform cell state annotation based on the feature gene list to obtain cell state annotation results to implement data classification for 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: performing filter selection 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 based on 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: Calculating the cell division rate based on the cell topology map and the pre-processed gene expression matrix to obtain cell division rate data; calculating 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 node weights according to the enhanced cell topology map to obtain initial node weights; Step S24: adjusting the initial node weights based on the cell division rate according to the enhanced cell topology map to obtain a division rate adjustment weight; adjusting the edge weights based on the RNA rate according to the enhanced cell topology map to obtain an RNA rate adjustment edge weight; 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 based on 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: Calculating the cell division rate based on 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: RNA rate inference is performed based on the mRNA matrix before and after splicing to obtain an RNA rate matrix; rate smoothing and noise reduction are performed 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 from 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 based on the node average RNA rate vector to obtain edge RNA rate difference data; Step S244: setting a minimum RNA rate difference according to the edge RNA rate difference data to obtain the minimum RNA rate difference data; Step S245: calculating the edge directional weight adjustment factor of the enhanced cell topology map according to the edge RNA rate difference data to obtain the directional weight adjustment factor; Step S246: performing edge weight adjustment on the edge RNA rate difference data according to the minimum RNA rate difference data and the directional weight adjustment factor to obtain an 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: performing fuzzy C-means clustering according to the trajectory feature matrix and the cluster quantity data to obtain the fuzzy clustering result; Step S34: constructing a membership matrix based on 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: performing clustering on 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 based on the gene expression grouping list and performing gene sorting to obtain gene entropy sorting data; Step S43: selecting a co-expression measurement method based on the pre-processed gene expression matrix to obtain a co-expression measurement method and parameters; Step S44: Calculating 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; performing network module identification on the gene co-expression network to obtain a network module; Step S46: Calculating the module information entropy of the network module according to the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy ranking data; 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.
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 based on the gene expression group list to obtain a kernel function; determining a bandwidth based on the gene expression group 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 bandwidth parameters to obtain gene cluster probability density; performing probability density normalization on the gene cluster probability density to obtain gene cluster probability distribution data; Step S423: Calculate gene information entropy based on gene clustering probability distribution data to obtain gene information entropy; Step S424: adjusting the gene information entropy weight 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 based on 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: Based on 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 based on the cell state prediction result to obtain a cell state annotation result.
10. A data classification system for single-cell sequencing, characterized in that: A data classification system for performing the single-cell sequencing data classification method according to claim 1, wherein the single-cell sequencing data classification system 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 based on the preprocessed gene expression matrix to obtain a cell topology map. A dynamic trajectory inference module is used to calculate the cell division rate and RNA rate based on the cell topology map and the preprocessed gene expression matrix to obtain RNA rate data; add node rates to the cell topology map based on the cell division rate data and the RNA rate data to obtain an enhanced cell topology map; adjust node weights and edge weights based on the enhanced cell topology map to obtain division rate adjustment weights and RNA rate adjustment edge weights; perform dynamic trajectory inference based on the enhanced cell topology map, division rate adjustment weights, and 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 the membership matrix to obtain the cell membership matrix; The characteristic gene screening module is used to calculate the clustering probability distribution based on 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 based on the cell membership matrix and the preprocessed gene expression matrix to obtain module information entropy sorting data; and screen characteristic genes based on 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 cell states based on the feature gene list and obtain cell state annotation results to implement data classification operations for single-cell sequencing.
Citation Information
Patent Citations
Web service hybrid evolutionary clustering method based on membrane computing
CN110659363A
Single cell trajectory inference method based on manifold learning and main curve
CN112768001A