Single-cell RNA sequencing data feature selection method and system based on deletion pattern

Abstract

The present invention discloses a method and system for feature selection of single-cell RNA sequencing data based on deletion patterns. This method targets the "dropout" phenomenon that is prevalent in scRNA-seq sequencing data by binarizing scRNA-seq sequencing data to represent the deletion pattern. Based on this pattern, the present invention proposes a feature selection algorithm for the fusion of mutually exclusive pairs of gene co-occurrence graphs. Compared with traditional scRNA-seq data analysis methods that solve the dropout problem by sampling, dimensionality reduction, etc., the present invention makes full use of the information in the dropout pattern to improve the accuracy and reliability of feature selection. This method provides strong technical support for the study of single-cell transcriptomics and has broad application prospects in multiple fields such as biology, medicine, and pharmacy, mainly including: understanding cellular heterogeneity, disease mechanisms and diagnosis, drug development, and personalized medicine.

Description

Technical Field

[0001] The present invention belongs to the field of single-cell RNA sequencing, and in particular relates to a method and system for single-cell RNA sequencing data feature selection based on deletion patterns. Background Art

[0002] Single-cell RNA sequencing (scRNA-seq) is a powerful technology that can measure gene expression at the single-cell level, revealing cell-to-cell heterogeneity. In bulk RNA sequencing of cell populations, important subpopulation structures and progression trajectories may be obscured due to the mixing of cell populations. However, scRNA-seq can capture these subtle differences, providing us with deeper biological insights.

[0003] Through scRNA-seq, researchers can understand the changes in gene expression of cells during development, differentiation, disease, etc., thereby revealing the subpopulation structure and progression trajectory of cells. This information is crucial for understanding biological processes, discovering new biomarkers, and developing new treatments. However, scRNA-seq experiments typically generate large amounts of data, containing genome-wide gene expression measurements of thousands or more individual cells, which poses challenges to the computational analysis and interpretation of the data. There are several reasons why the computational analysis of scRNA-seq data is challenging, such as high dimensionality, measurement noise, missing values, and size imbalance between rare and abundant populations. One of the most important challenges is the "dropout" phenomenon, which is an important feature in single-cell RNA sequencing (scRNA-seq). It refers to the situation where certain genes are expressed at low or moderate levels in one cell, but are not detected in another cell of the same type. This phenomenon is mainly caused by the following reasons: the total amount of mRNA in a single cell is very small, which may cause some genes with low expression to be missed during the RNA capture process; the imperfect mRNA capture efficiency, even if mRNA exists, it may not be successfully sequenced due to capture efficiency issues; the randomness of mRNA expression, even in the same type of cells, the expression of some genes may fluctuate. The dropout phenomenon will cause high sparsity of scRNA-seq data, that is, the proportion of zero values (i.e., gene expression is zero) in the data is very high. This sparsity not only increases the difficulty of data analysis, but may also obscure the true biological information. Therefore, when analyzing scRNA-seq data, special attention should be paid to the dropout phenomenon.

[0004] To address the dropout problem, researchers have developed a series of methods, including data preprocessing techniques (such as standardization, normalization, missing value interpolation, etc.), dimensionality reduction techniques (such as principal component analysis, t-SNE, etc.), clustering algorithms (such as K-means, hierarchical clustering, DBSCAN, etc.), and visualization tools. These methods mainly start from three aspects: using appropriate statistical models to estimate the true gene expression level to correct the underestimation caused by dropout; reducing the occurrence of dropout by increasing sequencing depth or improving RNA capture technology; using some specific algorithms or tools to identify and correct dropout events during the data analysis process, such as using machine learning-based methods to predict and fill missing values.

[0005] However, these methods primarily focus on how to better display and analyze scRNA-seq data, with less attention paid to how to effectively select biologically meaningful features from this data. Therefore, it is crucial to develop a method that can fully account for the dropout phenomenon and effectively select biologically meaningful features from scRNA-seq data. This approach can not only help researchers more accurately identify cell types, cell states, and key genes in biological processes, but also provide valuable clues for subsequent gene function analysis and disease prediction. Summary of the Invention

[0006] This paper proposes a feature selection method and system for single-cell RNA sequencing (scRNA-seq) data. Compared with traditional scRNA-seq data analysis methods (such as sampling and dimensionality reduction), this paper fully utilizes the information in the dropout pattern and makes an assumption: although the sparsity of scRNA-seq data is mainly caused by deletions caused by noise and randomness at the single-cell level, genes in the same pathway often show similar deletion patterns (i.e., binary zero / non-zero patterns) in different cell types. Based on this assumption, the binary information of scRNA-seq sequencing data is used to make subsequent analysis more straightforward. Then, based on the binary data, a gene co-occurrence map is constructed. The gene co-occurrence map describes the co-expression of different genes in a cell. Each node represents a gene, and the edges between nodes indicate that the two genes are detected to be expressed simultaneously under specific conditions or states. Such a co-occurrence map provides us with an intuitive display of the gene network structure, helping us to identify synergistic or mutually exclusive relationships between genes. Based on the gene co-occurrence map, a feature selection algorithm was invented that takes into account the information in the dropout pattern. Traditional feature selection methods often ignore this dropout phenomenon, but the present invention can effectively handle this situation and pay special attention to those gene pairs that show mutual exclusivity in the co-occurrence map, and integrate these mutually exclusive pairs into the feature selection process. Through this method, gene features that have an important impact on the research objectives can be selected more accurately. The accuracy and reliability of feature selection are improved. This improved feature selection method is of great significance for single-cell transcriptomics research. By applying this improved feature selection method, researchers can more accurately identify gene feature subsets associated with specific cell types or functions. This will help to better understand cellular heterogeneity, reveal the occurrence and development mechanisms of diseases, guide drug development and personalized medicine, and other issues.

[0007] In order to achieve the above technical objectives, the technical solutions of the present invention are as follows:

[0008] A method for feature selection of single-cell RNA sequencing data based on deletion patterns, comprising the following steps:

[0009] A. Construct a data preprocessing and binarization module to preprocess the raw scRNA-seq sequencing data and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout model;

[0010] B. Calculate the co-occurrence of genes in cells using a binary expression matrix based on a dropout pattern and the chi-square statistic, determine the gene co-occurrence matrix, construct a gene co-occurrence map, and calculate the gene mutual exclusion matrix using the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout pattern, determine mutually exclusive gene pairs, and obtain a gene co-occurrence weighted map;

[0011] C. Based on the first-order and second-order similarities of the gene co-occurrence graph nodes, the large-scale information network embedding model LINE is trained using the local information of the gene co-occurrence graph to generate graph node embedding vectors based on the local information of the gene co-occurrence weighted graph;

[0012] D. Walk the gene co-occurrence weighted graph based on the improved random walk algorithm to obtain the probability of nodes being visited in the gene co-occurrence weighted graph. Use the mutually exclusive gene pairs as prior knowledge to obtain the importance of the graph nodes based on the global information of the gene co-occurrence weighted graph.

[0013] E. Based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the average distance between genes in the embedding space is calculated using the Euclidean distance, and the preliminary feature gene selection set is determined by the established threshold;

[0014] F. Sort the preliminary feature gene selection set based on the node importance of the global information of the gene co-occurrence weighted graph, sequentially extract feature subsets of different sizes, use each feature subset to cluster all cells, and select the feature subset with the best clustering evaluation result as the final feature selection result.

[0015] Furthermore, the step A includes the following specific steps:

[0016] 1.1) Data preprocessing: Gene expression matrix X of size H×U raw As input, H and U are the number of rows and columns of the gene expression matrix, rows represent genes, and columns represent cells. Genes expressed in less than 6% of cells and genes expressed in more than 94% of cells are filtered out. The filtered gene expression matrix is X flt ;

[0017] 1.2) X flt Binarization is performed, and all zero elements in the binarized matrix are retained, and all non-zero elements are converted to 1, to obtain the binary expression matrix X based on the dropout mode N×L , N and L are the number of rows and columns of the binary expression matrix based on the dropout model, respectively.

[0018] Furthermore, the step B includes the following specific steps:

[0019] 2.1) Use the chi-square statistic to calculate the co-occurrence between each pair of genes. Based on the binary data of two genes g1 and g2, the chi-square value is calculated as

[0020]

[0021] A represents the number of cells in which both genes g1 and g2 were detected, B represents the number of cells in which g1 was detected but g2 was not detected, C represents the number of cells in which g1 was not detected but g2 was detected, and D represents the number of cells in which neither genes g1 nor g2 were detected.

[0022] 2.2) Determine the gene co-occurrence threshold, randomly arrange all genes in rows to generate ten random data sets, and record the highest chi-square value in the random data For the actual data before random permutation, the Chi-square value, calculate the chi-square value threshold threshold co =μ co +σ co , where μ co and σ co The actual data are lower than The mean and standard deviation of the chi-square value are used to determine the gene pairs with chi-square values higher than the chi-square value threshold as gene co-occurrence pairs, and the gene co-occurrence unweighted graph and edge set E are constructed;

[0023] 2.3) Calculate the weight of the node edge in the gene co-occurrence graph by the Sorensen-Dice coefficient to obtain the weight matrix W N×N , N is the number of rows of the binary expression matrix based on the dropout mode, and the gene co-occurrence weighted graph is constructed, W N×N The element W(i,j) in is the weight of the edge between two directly connected nodes i and j in the gene co-occurrence weighted graph. W(i,j) is calculated as follows:

[0024]

[0025] Where Γ(·) is the set of neighbor nodes of a node in the gene co-occurrence graph determined by the edge set E, |Γ(·)| represents the number of neighbor nodes, node i and node j are directly connected, and W(i, j) represents the weight of the edge between two directly connected nodes i and j;

[0026] 2.4) Calculate the mutual exclusion matrix M through matrix operations N×N :

[0027] M N×N =(1-X N×L )·X N×L T

[0028] Mutual exclusion matrix M N×N Each element M[i][j] in the matrix X represents the mutual exclusivity between gene i and gene j. N×L T is the binary gene expression matrix X N×L Transpose, calculate the mutual exclusion matrix M N×N The average value of all elements in μ mut and standard deviation σ mut , gene mutual exclusivity threshold threshold mut Established as:

[0029] threshold mut =μ mut +λ×σ mut

[0030] In the formula, λ represents the proportional coefficient, and its initial value is set to 1;

[0031] 2.5) Using threshold mut Binarize the mutual exclusion matrix and transform the mutual exclusion matrix M N×N The element value in is higher than the threshold mut The value of is set to 1, and the rest are set to 0.

[0032] In real datasets, completely mutually exclusive gene pairs are very rare, because gene expression is often affected by multiple factors, including environment, cell type, disease state, etc. Therefore, a more flexible and dynamic method is used to define and quantify the mutual exclusivity between genes, using the threshold mut Binarize the mutual exclusion matrix and transform the mutual exclusion matrix M N×N The element value in is higher than the threshold mut The value of is set to 1, and the rest are set to 0.

[0033] Furthermore, if the initial mutual exclusion matrix is asymmetric, the λ value is gradually increased with a step size of 0.01, and the mutual exclusion matrix is re-binarized. In each step, it is checked whether the new mutual exclusion matrix is symmetric. If not, the threshold is continued to be increased. If so, the iteration is stopped and the threshold at this time is regarded as the final threshold to obtain the final mutual exclusion matrix M N×N , matrix M N×N The gene pairs corresponding to the elements with a median value of 1 are determined to be mutually exclusive gene pairs.

[0034] Ideally, the mutual exclusion matrix should be symmetrical, because if gene A is mutually exclusive with gene B, then gene B should also be mutually exclusive with gene A. However, due to the influence of data noise and other factors, the actual mutual exclusion matrix may not be completely symmetrical, so an iterative method is used to update the threshold.

[0035] Furthermore, the step C includes the following specific steps:

[0036] 3.1) The first-order similarity of the model LINE is predefined;

[0037] First, the first-order similarity is defined as the proximity between directly connected nodes i and j in the gene co-occurrence weighted graph, that is, the weight W(i,j) of the edge in the gene co-occurrence weighted graph;

[0038] Secondly, initialize the d-dimensional embedding vector of each node in the gene co-occurrence weighted graph u1,u2,…,u d is the corresponding 1, 2...d-dimensional value in the node d-dimensional vector;

[0039] Next, based on the defined first-order similarity, a joint probability distribution is constructed between directly connected nodes i and j in the gene co-occurrence weighted graph, so that the probability value in the joint probability distribution is positively correlated with the first-order proximity between the nodes;

[0040] For each undirected edge (i, j) in the gene co-occurrence weighted graph, the joint probability distribution between nodes i and j is definition

[0041] for:

[0042]

[0043] Where, are the d-dimensional vector representations of node i and node j respectively, and Represented as the corresponding 1, 2...d-dimensional values in the d-dimensional vector of node i and node j, respectively, defining the empirical probability distribution for:

[0044]

[0045] Where W is the normalization factor, W = ∑ i,j∈E W(i,j);

[0046] 3.2) Model LINE second-order similarity is predefined;

[0047] First, the second-order similarity is defined as the similarity between two nodes that are not directly connected in the gene co-occurrence weighted graph through sharing neighbor nodes;

[0048] Secondly, initialize each node in the gene co-occurrence weighted graph as the d-dimensional embedding vector of its neighbor nodes

[0049] Next, based on the defined second-order similarity of nodes, the conditional probability distribution between non-adjacent nodes i and p in the gene co-occurrence weighted graph is constructed so that the conditional probability distribution can reflect the degree of second-order similarity between the node pairs;

[0050] Conditional probability distribution of non-adjacent nodes i and p in the gene co-occurrence weighted graph It means that when a given node i, the probability of its neighbor node j appearing is:

[0051]

[0052] Where V is the set of all nodes in the graph, is the d-dimensional vector representing node i itself, is the transpose of the d-dimensional vector representing node p as a neighbor node, It is represented as a d-dimensional vector when node k is a neighbor node, and the empirical probability distribution is defined as for:

[0053]

[0054] Where W(i,p) represents the weight of the edge between node i and node p, d i is the out-degree of node i, d i =∑ k∈Γ(i) W(i,k), W(i,k) represents the weight of the edge between node i and node k;

[0055] 3.3) Construct loss function;

[0056] Use KL divergence as a metric to minimize the empirical probability distribution of first-order similarity and the joint probability distribution The difference between them, the optimization loss function of the first-order similarity is defined as:

[0057]

[0058] Use KL divergence as a metric to minimize the empirical probability distribution of second-order similarity and conditional probability distribution The difference between them, the optimization loss function of the second-order similarity is defined as:

[0059]

[0060] 3.4) Using negative sampling and asynchronous stochastic gradient descent algorithm to optimize the second-order similarity loss function;

[0061] For each node in the gene co-occurrence weighted graph, determine the nodes directly connected to the node in the gene co-occurrence weighted graph as positive samples, and randomly select R nodes that are not directly connected to the node as negative samples. Use negative samples to optimize the loss function of the second-order similarity;

[0062] This method approximates the overall loss function by randomly extracting negative samples, avoiding the calculation of all non-existent edges, thereby significantly reducing the computational cost;

[0063] The loss function of the optimized second-order similarity is:

[0064]

[0065] Where, is the sigmoid function, x represents the variable of the sigmoid function, P(v) is the distribution of negative samples, v is the variable of the negative sample distribution, E r～P(v) [·] represents the expectation of node r sampled from the negative sample distribution;

[0066] 3.5) Model LINE training;

[0067] Set the total loss function O = O1 + O2, use the ASGD algorithm to iteratively update O, randomly select a positive sample and K negative samples to calculate the gradient of the total loss function, and use the following formula to update the d-dimensional embedding vector of node i in the next iteration based on the calculated gradient at the tth iteration:

[0068]

[0069] Where α is the learning rate, u i t and u i t+1 Represent the d-dimensional embedding vector of node i at the t-th and t+1-th iterations respectively; the embedding vector of each node in the gene co-occurrence weighted graph is continuously updated in each iteration until the maximum number of iterations is reached, and the node embedding vector based on the local information of the gene co-occurrence weighted graph in the low-dimensional space is iteratively updated.

[0070] Furthermore, the step D includes the following specific steps:

[0071] 4.1) Random walk gene co-occurrence weighted graph;

[0072] Randomly select a node in the gene co-occurrence weighted graph as the initial node start_node, set the number of random walk rounds J and the maximum number of steps per random walk num_step, initialize an empty path array path to record the movement path of the random walk, initialize an empty cell array exclusive_path to store the mutually exclusive pairs involved in the random walk process, set the weight decay factor decay_factor and initialize the decay matrix decay_matrix = M N×N , initialize the probability matrix probabilities:

[0073] probabilities=1-decay_matrix

[0074] Among them, M N×N represents the mutual exclusion matrix, N is the number of rows of the binary expression matrix based on the dropout mode;

[0075] 4.2) Get the neighbor node j of the current node, add the current node to the path array path, and record the current node's exclusive pair in the cell array exclusive_path. If the current node has no neighbor nodes, stop immediately and go to 4.3). If there are neighbor nodes, walk to the next node according to the corresponding probability in the probability matrix probabilities and go to 4.4).

[0076] 4.3) After updating the attenuation matrix and probability matrix, proceed to 4.4);

[0077] The influence of the current node's mutual exclusion on the information will gradually weaken as the number of random walk steps increases. For each step, the value of the decay matrix decay_matrix(i,j) is updated for each pair of nodes (i,j) in the cell array exclusive_path:

[0078] decay_matrix(i,j)=decay_matrix(i,j)×decay_factor

[0079] After updating the decay matrix decay_matrix, also update the probability matrix:

[0080] probabilities=1-decay_matrix

[0081] 4.4) If the number of random walk rounds and the maximum number of iteration steps have not been reached, return to 4.2). Otherwise, the random walk stops. The path array path generated by each random walk iteration is used to count the number of times all nodes are visited during the random walk and normalize the number. The normalized number of times is used as the node importance node_importances based on the global information of the gene co-occurrence weighted graph. Then, the nodes are sorted in descending order based on the size of the node importance node_importances to obtain sort.

[0082] Furthermore, the step E includes the following specific steps:

[0083] 5.1) In the embedding space, based on the node embedding vectors of the local information of the gene co-occurrence weighted graph, the distances between all genes are calculated to obtain the distance matrix D between genes, and then the average distance between genes and other genes is calculated;

[0084] Among them, the distance between node i and node j is:

[0085]

[0086] in, and are the d-dimensional embedding vectors of node i and node j based on the local information of the gene co-occurrence weighted graph after initialization, and They are and The value corresponding to the s-th dimension of the vector;

[0087] 5.2) Calculate the average μ0 and standard deviation σ0 of the average distances of all genes, then set the distance threshold to μ0-σ0. Based on the distance threshold, select genes whose average distances with other genes are greater than the distance threshold as elements of the feature selection set, forming a preliminary feature selection set δ.

[0088] Furthermore, the step F includes the following specific steps:

[0089] 6.1) Sort the genes in the preliminary feature set δ based on sort. After sorting is completed, extract feature subsets sequentially from the sorted list, sequentially extracting the first 5 elements, the first 10 elements, and the first 15 elements as feature subsets until all elements in the preliminary feature selection set are selected;

[0090] 6.2) Take each feature subset as the input feature set. For each input feature set, obtain all elements in the input feature set in the original gene expression matrix X. rawThe gene expression amount in is used as the gene expression feature vector of each cell in the original gene expression matrix, and K-means clustering is performed on all cells based on the gene expression feature vector of the cell;

[0091] 6.3) Use the mean square error (MSE) to determine the optimal number of clusters k0. MSE represents the mean square of the distance from each cell to the centroid of the cluster to which it belongs. Select different k1 values and calculate the MSE after clustering. Select the k1 value with the smallest MSE after clustering as the optimal number of clusters k0. The K-means clustering result corresponding to the optimal number of clusters k0 is the final clustering result.

[0092] The MSE calculation formula is as follows:

[0093]

[0094] Where n is the number of cell samples, C i is the cell set of the ith cluster, μ i is the centroid of the ith cluster, x1 is the gene element in the input feature set in the original gene expression matrix X raw The feature vector composed of the gene expression levels of all cells in ;

[0095] 6.4) Normalize the original gene expression matrix and perform linear regression with the final clustering results:

[0096]

[0097] In the formula is the gene expression feature vector of cell i predicted by the regression model, x i represents the final cluster label of cell i, β0 is the intercept term, β1 is the slope term, which represents the effect of cluster label on gene expression, ε i is the error term, y i is the true cell gene expression feature vector;

[0098] 6.5) Calculate the mean squared error (MSE) from the regression residuals regression :

[0099]

[0100] Where n is the number of cell samples;

[0101] This MSE regression It reflects the degree of fit between clustering results and gene expression;

[0102] 6.6) Calculate the MSE of the final clustering results for each different feature subset regression value, the minimum MSE regressionThe feature subset corresponding to the value is taken as the optimal feature set to obtain the final feature selection result.

[0103] Because it can achieve the most accurate clustering while retaining the most information.

[0104] Second, a single-cell RNA sequencing data feature selection system based on deletion patterns, comprising:

[0105] Data preprocessing and binarization module: preprocess the raw scRNA-seq sequencing data and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout model;

[0106] Gene co-occurrence weighted graph construction module: This module uses a binary expression matrix based on a dropout pattern and the chi-square statistic to calculate the co-occurrence of genes in cells, determines the gene co-occurrence matrix, constructs a gene co-occurrence graph, and calculates the gene mutual exclusion matrix through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout pattern, determines mutually exclusive gene pairs, and obtains a gene co-occurrence weighted graph;

[0107] Graph node embedding vector generation module: Based on the first-order and second-order similarities of gene co-occurrence graph nodes, the large-scale information network embedding model LINE is trained using the local information of the gene co-occurrence graph to generate graph node embedding vectors based on the local information of the gene co-occurrence weighted graph;

[0108] Graph node importance calculation module: This module walks the gene co-occurrence weighted graph based on an improved random walk algorithm to obtain the probability of a node being visited. It also uses mutually exclusive gene pairs as prior knowledge to obtain the graph node importance based on the global information of the gene co-occurrence weighted graph.

[0109] Preliminary feature gene selection module: Based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the Euclidean distance is used to calculate the average distance between genes in the embedding space, and the preliminary feature gene selection set is determined by the established threshold;

[0110] Feature subset optimal selection module: The preliminary feature gene selection set is sorted based on the node importance of the global information of the gene co-occurrence weighted graph, feature subsets of different sizes are extracted sequentially, all cells are clustered using each feature subset, and the feature subset with the best clustering evaluation result is selected as the final feature selection result.

[0111] In a third aspect, a computer-readable storage medium stores a computer program thereon, which, when executed by a central processing unit, implements the above-mentioned single-cell RNA sequencing data feature selection method based on deletion patterns.

[0112] Beneficial effects

[0113] The technical effect of the present invention is that it uses data preprocessing and binarization modules to remove noise and irrelevant information, thereby improving the accuracy of subsequent analysis. The binary expression matrix can highlight dropout information, facilitating subsequent feature selection and co-occurrence analysis. Using the chi-square statistic to calculate the co-occurrence of gene pairs can quantify the consistency of gene expression in different cells, providing strong evidence for revealing cell states. The determined gene co-occurrence matrix and gene mutual exclusion matrix can reveal the interactions and potential regulatory relationships between genes, providing a new perspective for understanding cell function and disease mechanisms. Using the local information of the gene co-occurrence map to train the LINE model, node embedding vectors reflecting the relationships between genes can be generated, providing a basis for subsequent feature selection and cluster analysis. The defined first-order similarity and second-order similarity can capture direct and indirect interactions between genes, improving the accuracy and interpretability of node embedding vectors. The improved random walk algorithm can fully utilize gene mutual exclusion pairs as prior knowledge to generate node feature importance based on global information, providing a more comprehensive basis for feature selection. Using the Euclidean distance to calculate the average distance between genes in the embedding space can quantify the similarities and differences between genes, providing a basis for preliminary feature selection. The threshold can be adjusted based on actual needs to balance the sensitivity and specificity of feature selection. By sorting and extracting feature subsets of different sizes and selecting the optimal feature subset based on clustering evaluation results, we can comprehensively evaluate the impact of different feature combinations on cell clustering, ensuring the effectiveness of the selected features in revealing cellular heterogeneity, disease mechanisms, and drug development, providing a scientific basis for final feature selection. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] Figure 1 This is a schematic diagram of the overall process of the single-cell RNA sequencing data feature selection method based on the dropout mode of the present invention;

[0115] Figure 2 The present invention generates a gene co-occurrence map and a gene mutually exclusive pair diagram based on the gene co-occurrence matrix and the mutual exclusion matrix; wherein (a) is the gene co-occurrence matrix and the gene co-occurrence map, and (b) is the gene mutual exclusion matrix and the gene mutually exclusive pair;

[0116] Figure 3 Schematic diagram of the LINE model and the first-order and second-order similarity definitions for the large network embedding model of the present invention; first-order proximity and second-order proximity are centered on node 4, first-order similarity represents the similarity between neighboring nodes 6 directly connected to node 4, and second-order similarity represents the similarity between node 4 and the neighboring nodes directly connected to node 6;

[0117] Figure 4Schematic diagram of the random walk iterative process with prior knowledge of mutually exclusive gene pairs;

[0118] Figure 5 Schematic diagram of the experimental comparison between the feature selection method of the present invention and other single-cell RNA sequencing data feature selection methods; (a), (b), and (c) are comparative schematic diagrams of the clustering evaluation indicators silhouette index, NMI, and ARI, respectively, and the horizontal axis in the figure represents different feature subsets. DETAILED DESCRIPTION

[0119] The method described in this embodiment is implemented using the MATLAB programming language and uses a publicly available human peripheral blood mononuclear cell (PBMC) dataset provided by 10X Genomics and sequenced using an Illumina NextSeq 500 sequencer. The dataset contains a total of 2,700 single cells and covers the expression information of 13,714 genes.

[0120] like Figure 1-Figure 5 As shown, an embodiment of the present invention provides a method for feature selection of single-cell RNA sequencing data based on deletion patterns, comprising the following steps:

[0121] A. Construct a data preprocessing and binarization module to preprocess the raw scRNA-seq sequencing data and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout model;

[0122] The step A comprises the following specific steps:

[0123] 1.1) Data preprocessing: Gene expression matrix X of size H×U raw As input, H and U are the number of rows and columns of the gene expression matrix, rows represent genes and columns represent cells, H = 13714, U = 2700, and genes expressed in less than 6% of cells and more than 94% of cells are filtered out. The gene expression matrix after filtering is X flt ;

[0124] 1.2) X flt Binarization is performed, and all zero elements in the binarized matrix are retained, and all non-zero elements are converted to 1, to obtain the binary expression matrix X based on the dropout mode N×L , N and L are the number of rows and columns of the binary expression matrix based on the dropout model, N = 13714, L = 2638.

[0125] B. Calculate the co-occurrence of genes in cells using a binary expression matrix based on a dropout pattern and the chi-square statistic, determine the gene co-occurrence matrix, construct a gene co-occurrence map, and calculate the gene mutual exclusion matrix using the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout pattern, determine mutually exclusive gene pairs, and obtain a gene co-occurrence weighted map;

[0126] Described step B comprises the following specific steps:

[0127] 2.1) Use the chi-square statistic to calculate the co-occurrence between each pair of genes. Based on the binary data of two genes g1 and g2, the chi-square value is calculated as

[0128]

[0129] A represents the number of cells in which both genes g1 and g2 were detected, B represents the number of cells in which g1 was detected but g2 was not detected, C represents the number of cells in which g1 was not detected but g2 was detected, and D represents the number of cells in which neither genes g1 nor g2 were detected.

[0130] 2.2) Determine the gene co-occurrence threshold, randomly arrange all genes in rows to generate ten random data sets, and record the highest chi-square value in the random data For the actual data before random permutation, the The chi-square value, mean μ co =0.1058 and standard deviation σ co =0.0321, calculate the chi-square value threshold co =μ co +σ co , that is, the chi-square value threshold co =0.1379, where μ co and σ co The actual data are lower than The mean and standard deviation of the chi-square value are used to determine the gene pairs with chi-square values higher than the chi-square value threshold as gene co-occurrence pairs, and the gene co-occurrence unweighted graph and edge set E are constructed;

[0131] 2.3) Calculate the weight of the node edge in the gene co-occurrence graph by the Sorensen-Dice coefficient to obtain the weight matrix W N×N , N is the number of rows of the binary expression matrix based on the dropout mode, and the gene co-occurrence weighted graph is constructed, W N×N The element W(i,j) in is the weight of the edge between two directly connected nodes i and j in the gene co-occurrence weighted graph. W(i,j) is calculated as follows:

[0132]

[0133] Where Γ(·) is the set of neighbor nodes of a node in the gene co-occurrence graph determined by the edge set E, |Γ(·)| represents the number of neighbor nodes, node i and node j are directly connected, and W(i, j) represents the weight of the edge between two directly connected nodes i and j;

[0134] 2.4) Calculate the mutual exclusion matrix M through matrix operations N×N :

[0135] M N×N =(1-X N×L )·X N×L T

[0136] Mutual exclusion matrix M N×N Each element M[i][j] in the matrix X represents the mutual exclusivity between gene i and gene j. N×L T is the binary gene expression matrix X N×L Transpose, calculate the mutual exclusion matrix M N×N The average value of all elements in μ mut =0.8113 and standard deviation σ mut =0.1797, gene mutual exclusivity threshold mut Established as:

[0137] threshold mut =μ mut +λ×σ mut

[0138] Where λ represents the proportional coefficient, and its initial value is set to 1. In real datasets, completely mutually exclusive gene pairs are very rare, because gene expression is usually affected by many factors, including environment, cell type, disease state, etc. Therefore, a more flexible and dynamic method is used to define and quantify the mutual exclusivity between genes, using the threshold mut Binarize the mutual exclusion matrix and transform the mutual exclusion matrix M N×N The element value in is higher than the threshold mutThe value of is set to 1, and the rest are set to 0. Ideally, the mutual exclusion matrix should be symmetrical, because if gene A is mutually exclusive to gene B, then gene B should also be mutually exclusive to gene A. However, due to the influence of data noise and other factors, the actual mutual exclusion matrix may not be completely symmetrical. Therefore, an iterative method is used to update the threshold. If the initial mutual exclusion matrix is not symmetrical, the λ value is gradually increased with a step size of 0.01, and the mutual exclusion matrix is re-binarized. In each step, check whether the new mutual exclusion matrix is symmetrical. If not, continue to increase the threshold. If so, stop the iteration, threshold mut =0.9910, and the threshold at this time is regarded as the final threshold, and the final mutual exclusion matrix M is obtained. N×N , matrix M N×N The gene pairs corresponding to the elements with a median value of 1 are determined to be mutually exclusive gene pairs. 2.5) Using the threshold mut Binarize the mutual exclusion matrix and transform the mutual exclusion matrix M N×N The element value in is higher than the threshold mut The value of is set to 1, and the rest are set to 0.

[0139] In real datasets, completely mutually exclusive gene pairs are very rare, because gene expression is often affected by multiple factors, including environment, cell type, disease state, etc. Therefore, a more flexible and dynamic method is used to define and quantify the mutual exclusivity between genes, using the threshold mut Binarize the mutual exclusion matrix and transform the mutual exclusion matrix M N×N The element value in is higher than the threshold mut The value of is set to 1, and the rest are set to 0.

[0140] Furthermore, if the initial mutual exclusion matrix is asymmetric, the λ value is gradually increased with a step size of 0.01, and the mutual exclusion matrix is re-binarized. In each step, it is checked whether the new mutual exclusion matrix is symmetric. If not, the threshold is continued to be increased. If so, the iteration is stopped and the threshold at this time is regarded as the final threshold to obtain the final mutual exclusion matrix M N×N , matrix M N×N The gene pairs corresponding to the elements with a median value of 1 are determined to be mutually exclusive gene pairs.

[0141] Ideally, the mutual exclusion matrix should be symmetrical, because if gene A is mutually exclusive with gene B, then gene B should also be mutually exclusive with gene A. However, due to the influence of data noise and other factors, the actual mutual exclusion matrix may not be completely symmetrical, so an iterative method is used to update the threshold.

[0142] C. Based on the first-order and second-order similarities of the gene co-occurrence graph nodes, the large-scale information network embedding model LINE is trained using the local information of the gene co-occurrence graph to generate graph node embedding vectors based on the local information of the gene co-occurrence weighted graph;

[0143] Described step C comprises the following specific steps:

[0144] 3.1) The first-order similarity of the model LINE is predefined;

[0145] First, the first-order similarity is defined as the proximity between directly connected nodes i and j in the gene co-occurrence weighted graph, that is, the weight W(i,j) of the edge in the gene co-occurrence weighted graph;

[0146] Secondly, initialize the d-dimensional embedding vector of each node in the gene co-occurrence weighted graph u1,u2,…,u d is the value of the corresponding 1, 2, ..., d dimension in the node d-dimensional vector; d = 128;

[0147] Next, based on the defined first-order similarity, a joint probability distribution is constructed between directly connected nodes i and j in the gene co-occurrence weighted graph, so that the probability value in the joint probability distribution is positively correlated with the first-order proximity between the nodes;

[0148] For each undirected edge (i, j) in the gene co-occurrence weighted graph, the joint probability distribution between nodes i and j is Defined as:

[0149]

[0150] Where, are the d-dimensional vector representations of node i and node j respectively, and Represented as the corresponding 1, 2...d-dimensional values in the d-dimensional vector of node i and node j, respectively, defining the empirical probability distribution for:

[0151]

[0152] Where W is the normalization factor, W = ∑ i,j∈E W(i,j);

[0153] 3.2) Model LINE second-order similarity is predefined;

[0154] First, the second-order similarity is defined as the similarity between two nodes that are not directly connected in the gene co-occurrence weighted graph through sharing neighbor nodes;

[0155] Secondly, initialize each node in the gene co-occurrence weighted graph as the d-dimensional embedding vector of its neighbor nodes

[0156] Next, based on the defined second-order similarity of nodes, the conditional probability distribution between non-adjacent nodes i and p in the gene co-occurrence weighted graph is constructed so that the conditional probability distribution can reflect the degree of second-order similarity between the node pairs;

[0157] Conditional probability distribution of non-adjacent nodes i and p in the gene co-occurrence weighted graph It means that when a given node i, the probability of its neighbor node j appearing is:

[0158]

[0159] Where V is the set of all nodes in the graph, is the d-dimensional vector representing node i itself, is the transpose of the d-dimensional vector representing node p as a neighbor node, It is represented as a d-dimensional vector when node k is a neighbor node, and the empirical probability distribution is defined as for:

[0160]

[0161] Where W(i,p) represents the weight of the edge between node i and node p, d i is the out-degree of node i, d i =∑ k∈Γ(i) W(i,k), W(i,k) represents the weight of the edge between node i and node k;

[0162] 3.3) Construct loss function;

[0163] Use KL divergence as a metric to minimize the empirical probability distribution of first-order similarity and the joint probability distribution The difference between them, the optimization loss function of the first-order similarity is defined as:

[0164]

[0165] Use KL divergence as a metric to minimize the empirical probability distribution of second-order similarity and conditional probability distribution The difference between them, the optimization loss function of the second-order similarity is defined as:

[0166]

[0167] 3.4) Using negative sampling and asynchronous stochastic gradient descent algorithm to optimize the second-order similarity loss function;

[0168] For each node in the gene co-occurrence weighted graph, determine the nodes directly connected to the node in the gene co-occurrence weighted graph as positive samples, and randomly select R nodes that are not directly connected to the node as negative samples. In this example, R = 50. Use negative samples to optimize the loss function of the second-order similarity;

[0169] This method approximates the overall loss function by randomly extracting negative samples, avoiding the calculation of all non-existent edges, thereby significantly reducing the computational cost;

[0170] The loss function of the optimized second-order similarity is:

[0171]

[0172] Where, is the sigmoid function, x represents the variable of the sigmoid function, P(v) is the distribution of negative samples, v is the variable of the negative sample distribution, E r～P(v) [·] represents the expectation of node r sampled from the negative sample distribution;

[0173] 3.5) Model LINE training;

[0174] Set the total loss function O = O1 + O2, use the ASGD algorithm to iteratively update O, randomly select a positive sample and K = 50 negative samples to calculate the gradient of the total loss function, and use the following formula to update the d-dimensional embedding vector of node i at the next iteration based on the calculated gradient at the tth iteration:

[0175]

[0176] Where α is the learning rate, u i t and u i t+1 Represent the d-dimensional embedding vector of node i at the t-th and t+1-th iterations respectively; the embedding vector of each node in the gene co-occurrence weighted graph is continuously updated in each iteration until the maximum number of iterations is reached, and the node embedding vector based on the local information of the gene co-occurrence weighted graph in the low-dimensional space is iteratively updated.

[0177] D. Walk the gene co-occurrence weighted graph based on the improved random walk algorithm to obtain the probability of nodes being visited in the gene co-occurrence weighted graph. Use the mutually exclusive gene pairs as prior knowledge to obtain the importance of the graph nodes based on the global information of the gene co-occurrence weighted graph.

[0178] Described step D comprises the following specific steps:

[0179] 4.1) Random walk gene co-occurrence weighted graph;

[0180] Randomly select a node in the gene co-occurrence weighted graph as the initial node start_node, set the number of random walk rounds J = 100 and the maximum number of steps per round of random walk num_step = 1000, initialize an empty path array path to record the movement path of the random walk, initialize an empty cell array exclusive_path to store the mutually exclusive pairs involved in the random walk process, set the weight decay factor decay_factor, in this example, according to the empirical value of 0.95 and initialize the decay matrix decay_matrix = M N×N , initialize the probability matrix probabilities:

[0181] probabilities=1-decay_matrix

[0182] Among them, M N×N represents the mutual exclusion matrix, N is the number of rows of the binary expression matrix based on the dropout mode;

[0183] 4.2) Get the neighbor node j of the current node, add the current node to the path array path, and record the current node's exclusive pair in the cell array exclusive_path. If the current node has no neighbor nodes, stop immediately and go to 4.3). If there are neighbor nodes, walk to the next node according to the corresponding probability in the probability matrix probabilities and go to 4.4).

[0184] 4.3) After updating the attenuation matrix and probability matrix, proceed to 4.4);

[0185] The influence of the current node's mutual exclusion on the information will gradually weaken as the number of random walk steps increases. For each step, the value of the decay matrix decay_matrix(i,j) is updated for each pair of nodes (i,j) in the cell array exclusive_path:

[0186] decay_matrix(i,j)=decay_matrix(i,j)×decay_factor

[0187] After updating the decay matrix decay_matrix, also update the probability matrix:

[0188] probabilities=1-decay_matrix

[0189] 4.4) If the number of random walk rounds and the maximum number of iteration steps have not been reached, return to 4.2). Otherwise, the random walk stops. The path array path generated by each random walk iteration is used to count the number of times all nodes are visited during the random walk and normalize the number. The normalized number of times is used as the node importance node_importances based on the global information of the gene co-occurrence weighted graph. Then, the nodes are sorted in descending order based on the size of the node importance node_importances to obtain sort.

[0190] E. Based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the average distance between genes in the embedding space is calculated using the Euclidean distance, and the preliminary feature gene selection set is determined by the established threshold;

[0191] The step E comprises the following specific steps:

[0192] 5.1) In the embedding space, based on the node embedding vectors of the local information of the gene co-occurrence weighted graph, the distances between all genes are calculated to obtain the distance matrix D between genes, and then the average distance between genes and other genes is calculated;

[0193] Among them, the distance between node i and node j is:

[0194]

[0195] in, and are the d-dimensional embedding vectors of node i and node j based on the local information of the gene co-occurrence weighted graph after initialization, and They are and The value corresponding to the s-th dimension of the vector;

[0196] 5.2) Calculate the average μ0 and standard deviation σ0 of the average distances of all genes, then set the distance threshold to μ0-σ0. Based on the distance threshold, select genes whose average distances with other genes are greater than the distance threshold as elements of the feature selection set, forming a preliminary feature selection set δ.

[0197] F. Sort the preliminary feature gene selection set based on the node importance of the global information of the gene co-occurrence weighted graph, sequentially extract feature subsets of different sizes, use each feature subset to cluster all cells, and select the feature subset with the best clustering evaluation result as the final feature selection result.

[0198] The step F comprises the following specific steps:

[0199] 6.1) Sort the genes in the preliminary feature set δ based on sort. After sorting is completed, extract feature subsets sequentially from the sorted list, sequentially extracting the first 5 elements, the first 10 elements, and the first 15 elements as feature subsets until all elements in the preliminary feature selection set are selected;

[0200] 6.2) Take each feature subset as the input feature set. For each input feature set, obtain all elements in the input feature set in the original gene expression matrix X. raw The gene expression values in the original gene expression matrix are used as the gene expression feature vectors of each cell in the original gene expression matrix. K-means clustering is performed on all cells based on the gene expression feature vectors of the cells. The goal of K-means clustering is to divide the 13714 cells into 12 clusters.

[0201] 6.3) Use the mean square error (MSE) to determine the optimal number of clusters k0. MSE represents the mean square of the distance from each cell to the centroid of the cluster to which it belongs. Select different k1 values and calculate the MSE after clustering. Select the k1 value with the smallest MSE after clustering as the optimal number of clusters k0. The K-means clustering result corresponding to the optimal number of clusters k0 is the final clustering result.

[0202] The MSE calculation formula is as follows:

[0203]

[0204] Where n is the number of cell samples, C i is the cell set of the ith cluster, μ i is the centroid of the ith cluster, x1 is the gene element in the input feature set in the original gene expression matrix X raw The feature vector composed of the gene expression levels of all cells in ;

[0205] 6.4) Normalize the original gene expression matrix and perform linear regression with the final clustering results:

[0206]

[0207] In the formula is the gene expression feature vector of cell i predicted by the regression model, x i represents the final cluster label of cell i, β0 is the intercept term, β1 is the slope term, which represents the effect of cluster label on gene expression, ε i is the error term, y i is the true cell gene expression feature vector;

[0208] 6.5) Calculate the mean squared error (MSE) from the regression residuals regression :

[0209]

[0210] Where n is the number of cell samples;

[0211] This MSE regression It reflects the degree of fit between clustering results and gene expression;

[0212] 6.6) Calculate the MSE of the final clustering results for each different feature subset regression value, the minimum MSE regression The feature subset corresponding to the value is taken as the optimal feature set to obtain the final feature selection result.

[0213] Because it can achieve the most accurate clustering while retaining the most information.

[0214] The single-cell RNA sequencing data feature selection method based on the dropout mode of the present invention is compared with three other feature selection methods based on highly variable genes. The three methods based on highly variable genes involved in the comparison are HVG, FEAST and DUBStepR. It can be seen from the clustering result comparison chart that: for the Silhouette Index indicator, as the size of the feature subset changes, the broken line of the method of the present invention is always at a higher position, and when the feature subset is within 150 genes, its value is significantly higher than the other three methods. This shows that the method of the present invention performs well in measuring the degree of closeness and separation between samples, and can better distinguish different cell types; for the NMI (Normalized Mutual Information) and ARI (Adjusted Rand Index) indicators, the method of the present invention also has good performance under each feature set size. Compared with the other three methods, it can more accurately reflect the consistency between the clustering results and the true label, especially when the feature set size is 100, the advantage is more obvious. This means that the method of the present invention has higher reliability in evaluating the accuracy and stability of clustering results.

[0215] The present invention also provides a single-cell RNA sequencing data feature selection system based on deletion patterns, comprising:

[0216] Data preprocessing and binarization module: preprocess the raw scRNA-seq sequencing data and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout model;

[0217] Gene co-occurrence weighted graph construction module: This module uses a binary expression matrix based on a dropout pattern and the chi-square statistic to calculate the co-occurrence of genes in cells, determines the gene co-occurrence matrix, constructs a gene co-occurrence graph, and calculates the gene mutual exclusion matrix through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout pattern, determines mutually exclusive gene pairs, and obtains a gene co-occurrence weighted graph;

[0218] Graph node embedding vector generation module: Based on the first-order and second-order similarities of gene co-occurrence graph nodes, the large-scale information network embedding model LINE is trained using the local information of the gene co-occurrence graph to generate graph node embedding vectors based on the local information of the gene co-occurrence weighted graph;

[0219] Graph node importance calculation module: This module walks the gene co-occurrence weighted graph based on an improved random walk algorithm to obtain the probability of a node being visited. It also uses mutually exclusive gene pairs as prior knowledge to obtain the graph node importance based on the global information of the gene co-occurrence weighted graph.

[0220] Preliminary feature gene selection module: Based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the Euclidean distance is used to calculate the average distance between genes in the embedding space, and the preliminary feature gene selection set is determined by the established threshold;

[0221] Feature subset optimal selection module: The preliminary feature gene selection set is sorted based on the node importance of the global information of the gene co-occurrence weighted graph, feature subsets of different sizes are extracted sequentially, all cells are clustered using each feature subset, and the feature subset with the best clustering evaluation result is selected as the final feature selection result.

[0222] It should be understood that the specific implementation process of the above-mentioned unit modules is referred to in the method content, and the present invention does not provide a detailed description here. Moreover, the division of the above-mentioned functional module units is merely a division of logical functions. In actual implementation, other division methods may be used. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not implemented. At the same time, the above-mentioned integrated units can be implemented in the form of hardware or software functional units.

[0223] The present invention also provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a central processing unit, the computer program implements the above-mentioned single-cell RNA sequencing data feature selection method based on deletion patterns.

[0224] For the specific implementation process, please refer to the description of the above method.

[0225] The readable storage medium is a computer-readable storage medium, which may be an internal storage unit of the controller described in any of the aforementioned embodiments, such as a hard disk or memory of the controller. The readable storage medium may also be an external storage device of the controller, such as a plug-in hard disk, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), etc. equipped on the controller. Furthermore, the readable storage medium may also include both an internal storage unit of the controller and an external storage device. The readable storage medium is used to store the computer program and other programs and data required by the controller. The readable storage medium may also be used to temporarily store data that has been output or is to be output.

[0226] Based on this understanding, the technical solution of the present invention, or the portion that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes instructions for causing a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present invention. The aforementioned readable storage medium includes various media capable of storing program code, such as a USB flash drive, a mobile hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk.

[0227] It should be emphasized that the examples described in the present invention are illustrative rather than restrictive. Therefore, the present invention is not limited to the examples described in the specific embodiments. Any other embodiments derived by those skilled in the art based on the technical solution of the present invention that do not depart from the purpose and scope of the present invention, whether modified or replaced, also fall within the scope of protection of the present invention.

Claims

1. A method for feature selection of single-cell RNA sequencing data based on deletion patterns, characterized in that: The following steps are involved: A. Construct a data preprocessing and binarization module to preprocess the raw scRNA-seq sequencing data and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout model; B. Calculate the co-occurrence of genes in cells using a binary expression matrix based on a dropout pattern and the chi-square statistic, determine the gene co-occurrence matrix, construct a gene co-occurrence map, and calculate the gene mutual exclusion matrix using the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout pattern, determine mutually exclusive gene pairs, and obtain a gene co-occurrence weighted map; C. Based on the first-order and second-order similarities of the gene co-occurrence graph nodes, the large-scale information network embedding model LINE is trained using the local information of the gene co-occurrence graph to generate graph node embedding vectors based on the local information of the gene co-occurrence weighted graph; D. Walk the gene co-occurrence weighted graph based on the improved random walk algorithm to obtain the probability of nodes being visited in the gene co-occurrence weighted graph. Use the mutually exclusive gene pairs as prior knowledge to obtain the importance of the graph nodes based on the global information of the gene co-occurrence weighted graph. E. Based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the average distance between genes in the embedding space is calculated using the Euclidean distance, and the preliminary feature gene selection set is determined by the established threshold; F. Sort the preliminary feature gene selection set based on the node importance of the global information of the gene co-occurrence weighted graph, sequentially extract feature subsets of different sizes, use each feature subset to cluster all cells, and select the feature subset with the best clustering evaluation result as the final feature selection result; Described step D comprises the following specific steps: 4.1) Random walk gene co-occurrence weighted graph; Randomly select a node in the gene co-occurrence weighted graph as the initial node start_node, set the number of random walk rounds J and the maximum number of steps per random walk num_step, initialize an empty path array path to record the movement path of the random walk, initialize an empty cell array exclusive_path to store the mutually exclusive pairs involved in the random walk process, set the weight decay factor decay_factor and initialize the decay matrix decay_matrix = M N×N , initialize the probability matrix probabilities: probabilities=1-decay_matrix; Among them, M N×N represents the mutual exclusion matrix, N is the number of rows of the binary expression matrix based on the dropout mode; 4.2) Get the neighbor node j of the current node, add the current node to the path array path, and record the current node's exclusive pair in the cell array exclusive_path. If the current node has no neighbor nodes, stop immediately and go to 4.3). If there are neighbor nodes, walk to the next node according to the corresponding probability in the probability matrix probabilities and go to 4.4). 4.3) After updating the attenuation matrix and probability matrix, proceed to 4.4); The influence of the current node's mutual exclusion on the information will gradually weaken as the number of random walk steps increases. For each step, the value of the decay matrix decay_matrix(i,j) is updated for each pair of nodes (i,j) in the cell array exclusive_path: decay_matrix(i,j)=decay_matrix(i,j)×decay_factor; After updating the decay matrix decay_matrix, also update the probability matrix: probabilities=1-decay_matrix; 4.4) If the number of random walk rounds and the maximum number of iteration steps have not been reached, return to 4.2). Otherwise, the random walk stops. The path array path generated by each random walk iteration is used to count the number of times all nodes are visited during the random walk and normalize the number. The normalized number of times is used as the node importance node_importances based on the global information of the gene co-occurrence weighted graph. Then, the nodes are sorted in descending order based on the size of the node importance node_importances to obtain sort.

2. The method according to claim 1, characterized in that The step A comprises the following specific steps: 1.1) Data preprocessing: Gene expression matrix X of size H×U raw As input, H and U are the number of rows and columns of the gene expression matrix, rows represent genes, and columns represent cells. Genes expressed in less than 6% of cells and genes expressed in more than 94% of cells are filtered out. The filtered gene expression matrix is X flt ; 1.2) X flt Binarization is performed, and all zero elements in the binarized matrix are retained, and all non-zero elements are converted to 1, to obtain the binary expression matrix X based on the dropout mode N×L , N and L are the number of rows and columns of the binary expression matrix based on the dropout model, respectively.

3. The method according to claim 1, characterized in that Described step B comprises the following specific steps: 2.1) Use the chi-square statistic to calculate the co-occurrence between each pair of genes. Based on the binary data of two genes g1 and g2, the chi-square value is calculated as: Among them, A represents the number of cells in which both genes g1 and g2 were detected, B represents the number of cells in which g1 was detected but g2 was not detected, C represents the number of cells in which g1 was not detected but g2 was detected, and D represents the number of cells in which neither genes g1 nor g2 were detected. 2.2) Determine the gene co-occurrence threshold, randomly arrange all genes in rows to generate ten random data sets, and record the highest chi-square value in the random data For the actual data before random permutation, the Chi-square value, calculate the chi-square value threshold threshold co =μ co +σ co , where μ co and σ co The actual data are lower than The mean and standard deviation of the chi-square value are used to determine the gene pairs with chi-square values higher than the chi-square value threshold as gene co-occurrence pairs, and the gene co-occurrence unweighted graph and edge set E are constructed; 2.3) Calculate the weight of the node edge in the gene co-occurrence graph by the Sorensen-Dice coefficient to obtain the weight matrix W N×N , N is the number of rows of the binary expression matrix based on the dropout mode, and the gene co-occurrence weighted graph is constructed, W N×N The element W(i,j) in is the weight of the edge between two directly connected nodes i and j in the gene co-occurrence weighted graph. W(i,j) is calculated as follows: Where Γ(·) is the set of neighbor nodes of a node in the gene co-occurrence graph determined by the edge set E, |Γ(·)| represents the number of neighbor nodes, and node i and node j are directly connected; 2.4) Calculate the mutual exclusion matrix M through matrix operations N×N : M N×N =(1-X N×L )·X N×L T Mutual Exclusion Matrix M N×N Each element M[i][j] in the matrix X represents the mutual exclusivity between gene i and gene j. N×L T is the binary gene expression matrix X N×L Transpose, calculate the mutual exclusion matrix M N×N The average value of all elements in μ mut and standard deviation σ mut , gene mutual exclusivity threshold threshold mut Established as: threshold mut =μ mut +λ×σ mut In the formula, λ represents the proportional coefficient, and its initial value is set to 1; 2.5) Using threshold mut Binarize the mutual exclusion matrix and transform the mutual exclusion matrix M N×N The element value in is higher than the threshold mut The value of is set to 1, and the rest are set to 0.

4. The method according to claim 3, characterized in that If the initial mutual exclusion matrix is asymmetric, gradually increase the λ value with a step size of 0.01 and re-binarize the mutual exclusion matrix. In each step, check whether the new mutual exclusion matrix is symmetric. If not, continue to increase the threshold. If so, stop the iteration and regard the threshold at this time as the final threshold to obtain the final mutual exclusion matrix M N×N , M N×N The gene pairs corresponding to the elements with a median value of 1 are determined to be mutually exclusive gene pairs.

5. The method according to claim 3, characterized in that Described step C comprises the following specific steps: 3.1) The first-order similarity of the model LINE is predefined; First, the first-order similarity is defined as the proximity between directly connected nodes i and j in the gene co-occurrence weighted graph, that is, the weight W(i,j) of the edge in the gene co-occurrence weighted graph; Secondly, initialize the d-dimensional embedding vector of each node in the gene co-occurrence weighted graph u1,u2,…,u d is the corresponding 1, 2...d-dimensional value in the node d-dimensional vector; Next, based on the defined first-order similarity, a joint probability distribution is constructed between directly connected nodes i and j in the gene co-occurrence weighted graph, so that the probability value in the joint probability distribution is positively correlated with the first-order proximity between the nodes; For each undirected edge (i, j) in the gene co-occurrence weighted graph, the joint probability distribution between nodes i and j is Defined as: Where, are the d-dimensional vector representations of node i and node j respectively, and Represented as the corresponding 1, 2...d-dimensional values in the d-dimensional vector of node i and node j, respectively, defining the empirical probability distribution for: Where W is the normalization factor, W = Σ i,j∈E W(i,j); 3.2) Model LINE second-order similarity is predefined; First, the second-order similarity is defined as the similarity between two nodes that are not directly connected in the gene co-occurrence weighted graph through sharing neighbor nodes; Secondly, initialize each node in the gene co-occurrence weighted graph as the d-dimensional embedding vector of its neighbor nodes Next, based on the defined second-order similarity of nodes, the conditional probability distribution between non-adjacent nodes i and p in the gene co-occurrence weighted graph is constructed so that the conditional probability distribution can reflect the degree of second-order similarity between the node pairs; Conditional probability distribution of non-adjacent nodes i and p in the gene co-occurrence weighted graph It means that when a given node i, the probability of its neighbor node j appearing is: Where V is the set of all nodes in the graph, is the d-dimensional vector representing node i itself, is the transpose of the d-dimensional vector representing node p as a neighbor node, It is represented as a d-dimensional vector when node k is a neighbor node, and the empirical probability distribution is defined as for: Where W(i,p) represents the weight of the edge between node i and node p, d i is the out-degree of node i, d i =Σ k∈Γ(i) W(i,k), W(i,k) represents the weight of the edge between node i and node k; 3.3) Construct loss function; Use KL divergence as a metric to minimize the empirical probability distribution of first-order similarity and the joint probability distribution The difference between them, the optimization loss function of the first-order similarity is defined as: Use KL divergence as a metric to minimize the empirical probability distribution of second-order similarity and conditional probability distribution The difference between them, the optimization loss function of the second-order similarity is defined as: 3.4) Using negative sampling and asynchronous stochastic gradient descent algorithm to optimize the second-order similarity loss function; For each node in the gene co-occurrence weighted graph, determine the nodes directly connected to the node in the gene co-occurrence weighted graph as positive samples, and randomly select R nodes that are not directly connected to the node as negative samples. Use negative samples to optimize the loss function of the second-order similarity; The loss function of the optimized second-order similarity is: Where, is the sigmoid function, x represents the variable of the sigmoid function, P(v) is the distribution of negative samples, v is the variable of the negative sample distribution, E r～P(v) [·] represents the expectation of node r sampled from the negative sample distribution; 3.5) Model LINE training; Set the total loss function O = O1 + O2, use the ASGD algorithm to iteratively update O, randomly select a positive sample and K negative samples to calculate the gradient of the total loss function, and use the following formula to update the d-dimensional embedding vector of node i in the next iteration based on the calculated gradient at the tth iteration: Where α is the learning rate, u i t and u i t+1 Represent the d-dimensional embedding vector of node i at the t-th and t+1-th iterations respectively; the embedding vector of each node in the gene co-occurrence weighted graph is continuously updated in each iteration until the maximum number of iterations is reached, and the node embedding vector based on the local information of the gene co-occurrence weighted graph in the low-dimensional space is iteratively updated.

6. The method according to claim 5, characterized in that The step E comprises the following specific steps: 5.1) In the embedding space, based on the node embedding vectors of the local information of the gene co-occurrence weighted graph, the distances between all genes are calculated to obtain the distance matrix D between genes, and then the average distance between genes and other genes is calculated; Among them, the distance between node i and node j is: in, and are the d-dimensional embedding vectors of node i and node j based on the local information of the gene co-occurrence weighted graph after initialization, and They are and The value corresponding to the s-th dimension of the vector; 5.2) Calculate the average μ0 and standard deviation σ0 of the average distances of all genes, then set the distance threshold to μ0-σ0. Based on the distance threshold, select genes whose average distances with other genes are greater than the distance threshold as elements of the feature selection set, forming a preliminary feature selection set δ.

7. The method according to claim 1, characterized in that The step F comprises the following specific steps: 6.1) Sort the genes in the preliminary feature set δ based on sort. After sorting is completed, extract feature subsets sequentially from the sorted list, sequentially extracting the first 5 elements, the first 10 elements, and the first 15 elements as feature subsets until all elements in the preliminary feature selection set are selected; 6.2) Take each feature subset as the input feature set. For each input feature set, obtain all elements in the input feature set in the original gene expression matrix X. raw The gene expression amount in is used as the gene expression feature vector of each cell in the original gene expression matrix, and K-means clustering is performed on all cells based on the gene expression feature vector of the cell; 6.3) Use the mean square error (MSE) to determine the optimal number of clusters k0. MSE represents the mean square of the distance from each cell to the centroid of the cluster to which it belongs. Select different k1 values and calculate the MSE after clustering. Select the k1 value with the smallest MSE after clustering as the optimal number of clusters k0. The K-means clustering result corresponding to the optimal number of clusters k0 is the final clustering result. The MSE calculation formula is as follows: Where n is the number of cell samples, C i is the cell set of the ith cluster, μ i is the centroid of the ith cluster, x1 is the gene element in the input feature set in the original gene expression matrix X raw The feature vector composed of the gene expression levels of all cells in ; 6.4) Normalize the original gene expression matrix and perform linear regression with the final clustering results: In the formula is the gene expression feature vector of cell i predicted by the regression model, x i represents the final cluster label of cell i, β0 is the intercept term, β1 is the slope term, which represents the effect of cluster label on gene expression, ε i is the error term, y i is the true cell gene expression feature vector; 6.5) Calculate the mean squared error (MSE) from the regression residuals regression : Where n is the number of cell samples; 6.6) Calculate the MSE of the final clustering results for each different feature subset regression value, the minimum MSE regression The feature subset corresponding to the value is taken as the optimal feature set to obtain the final feature selection result.

8. A single-cell RNA sequencing data feature selection system based on deletion patterns, characterized in that: include: Data preprocessing and binarization module: preprocess the raw scRNA-seq sequencing data and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout model; Gene co-occurrence weighted graph construction module: This module uses a binary expression matrix based on a dropout pattern and the chi-square statistic to calculate the co-occurrence of genes in cells, determines the gene co-occurrence matrix, constructs a gene co-occurrence graph, and calculates the gene mutual exclusion matrix through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout pattern, determines mutually exclusive gene pairs, and obtains a gene co-occurrence weighted graph; Graph node embedding vector generation module: Based on the first-order and second-order similarities of gene co-occurrence graph nodes, the large-scale information network embedding model LINE is trained using the local information of the gene co-occurrence graph to generate graph node embedding vectors based on the local information of the gene co-occurrence weighted graph; Graph node importance calculation module: This module walks the gene co-occurrence weighted graph based on an improved random walk algorithm to obtain the probability of a node being visited. It also uses mutually exclusive gene pairs as prior knowledge to obtain the graph node importance based on the global information of the gene co-occurrence weighted graph. The graph node importance calculation module uses the single-cell RNA sequencing data feature selection method based on deletion patterns according to claim 1 to calculate the graph node importance; Preliminary feature gene selection module: Based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the Euclidean distance is used to calculate the average distance between genes in the embedding space, and the preliminary feature gene selection set is determined by the established threshold; Feature subset optimal selection module: The preliminary feature gene selection set is sorted based on the node importance of the global information of the gene co-occurrence weighted graph, feature subsets of different sizes are extracted sequentially, all cells are clustered using each feature subset, and the feature subset with the best clustering evaluation result is selected as the final feature selection result.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When executed by a central processing unit, the computer program implements a single-cell RNA sequencing data feature selection method based on deletion patterns as described in any one of claims 1 to 7.

Images