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

Abstract

The invention discloses a single-cell RNA (Ribonucleic Acid) sequencing data feature selection method and a single-cell RNA sequencing data feature selection system based on a deletion pattern, and the method is used for representing the deletion pattern by binarizing scRNA-seq sequencing data aiming at a'dropout 'phenomenon generally existing in the scRNA-seq sequencing data. Based on the mode, the invention provides a feature selection algorithm for gene co-occurrence graph mutual exclusion pair fusion. Compared with a traditional scRNA-seq data analysis method for solving the falling problem through sampling, dimensionality reduction and the like, information in a dropout mode is fully utilized, and the accuracy and reliability of feature selection are improved. The method provides powerful technical support for research of single cell transcriptomics, and has wide application prospects in multiple fields of biology, medicine, pharmacy and the like, mainly including understanding of cell heterogeneity, disease mechanism and diagnosis, drug research and development, personalized medical treatment and the like.

Description

Technical Field

[0001] The present invention belongs to the field of single-cell RNA sequencing, and in particular relates to a single-cell RNA sequencing data feature selection method and system 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, thereby revealing cell-to-cell heterogeneity. In bulk cell population RNA sequencing, some important subpopulation structures and progression trajectories may be obscured due to the mixing of cell populations. However, scRNA-seq can capture these subtle differences and provide 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 usually generate a large amount 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), which 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 there is mRNA, it may not be successfully sequenced due to capture efficiency problems; 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., zero gene expression) in the data is very high. This sparsity not only increases the difficulty of data analysis, but may also obscure the real biological information. Therefore, when analyzing scRNA-seq data, special attention should be paid to the dropout phenomenon.

[0004] In order to deal with 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, etc. 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; in the process of data analysis, using some specific algorithms or tools to identify and correct dropout events, such as using machine learning-based methods to predict and fill missing values.

[0005] However, these methods mainly focus on how to better display and analyze scRNA-seq data, but pay less attention to how to effectively select biologically meaningful features from these data. Therefore, it is of great significance to develop a method that can fully consider the dropout phenomenon and effectively select biologically meaningful features from scRNA-seq data. This method 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] The present invention 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, dimensionality reduction, etc.), the present invention makes full use of the information in the dropout pattern and makes an assumption: although the sparsity of scRNA-seq data is mainly caused by the loss caused by noise and randomness at the single-cell level, genes on the same pathway often show similar loss 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 direct. 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 cells, each node represents a gene, and the edges between the 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, which helps us identify the synergy or mutual exclusion relationship between genes. Based on the gene co-occurrence map, a feature selection algorithm was invented, which takes into account the information in the dropout pattern. Traditional feature selection methods often ignore this dropout phenomenon, while 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 the study of single-cell transcriptomics. 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 cell 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 solution of the present invention is as follows:

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

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

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

[0011] C. Based on the first-order similarity and second-order similarity 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 the gene co-occurrence weighted graph node being visited, and use the gene mutually exclusive pairs as prior knowledge to obtain the importance of the graph node based on the global information of the gene co-occurrence weighted graph;

[0013] E. The node embedding vector based on the local information of the gene co-occurrence weighted graph uses the Euclidean distance to calculate the average distance between genes in the embedding space, 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, extract feature subsets of different sizes in sequence, 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 comprises the following specific steps:

[0016] 1.1) Data preprocessing: Gene expression matrix X of size H×U raw is the input, H and U are the number of rows and columns of the gene expression matrix, rows represent genes, columns represent cells, and genes expressed in less than 6% of cells and in more than 94% of cells are filtered out. The gene expression matrix after filtering 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 mode, respectively.

[0018] Furthermore, the step B comprises 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 the 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 weights of the node edges in the gene co-occurrence graph using 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 a 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 The transpose of , calculate the mutual exclusion matrix M N×N The average value of all elements in μ mut and standard deviation σ mut , gene mutual exclusivity threshold mut Established as:

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

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

[0031] 2.5) Using threshold mut Binarize the mutual exclusion matrix and convert 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 usually 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 convert 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 yes, 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 , the matrix M N×N The gene pairs corresponding to the elements with a median value of 1 are determined as mutually exclusive gene pairs.

[0034] 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 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 comprises 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 u 1 ,u 2 ,…,u d is the corresponding 1, 2...d-dimensional value in the node's d-dimensional vector;

[0039] Next, based on the defined first-order similarity, the joint probability distribution between directly connected nodes i and j in the gene co-occurrence weighted graph is constructed, 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] In the formula, 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, define the empirical probability distribution for:

[0044]

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

[0046] 3.2) The second-order similarity of the model LINE 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 shared neighbor nodes;

[0048] Secondly, initialize each node in the gene co-occurrence weighted graph as the d-dimensional embedding vector of its neighboring 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 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 defines the optimization loss function of the first-order similarity 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 defines the optimization loss function of the second-order similarity 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, randomly select R nodes that are not directly connected to the node as negative samples, and use negative samples to optimize the loss function of the second-order similarity;

[0062] This method obtains negative samples by random sampling to approximate the overall loss function, avoiding the calculation of all non-existent edges, thus significantly reducing the computational cost;

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

[0064]

[0065] In the formula, 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 = O 1 +O 2 , 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] In the formula, α is the learning rate, u i t and u i t+1 They 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 obtained by iterative update.

[0070] Furthermore, the step D comprises 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 in each round of random walk num_step, initialize an empty path array path to record the moving 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 mutually exclusive 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 exclusive pair of the current node in the cell array exclusive_path. If the current node has no neighbor node, stop directly and go to 4.3). If there is a neighbor node, 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 mutual exclusion of the current node on the information will gradually weaken as the number of steps of the random walk increases. For each step, for each pair of nodes (i, j) in the cell array exclusive_path, the value of the decay matrix decay_matrix(i, j) is updated:

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

[0079] After updating the decay matrix decay_matrix, the probability matrix is ​​also updated:

[0080] probabilities=1-decay_matrix

[0081] 4.4) If the number of random walk rounds and the maximum number of iteration steps are not reached, return to 4.2), otherwise, the random walk stops, and 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 them. 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, and then the nodes are arranged in descending order based on the size of the node importance node_importances to obtain the sort.

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

[0083] 5.1) In the embedding space, based on the node embedding vector of the local information of the gene co-occurrence weighted graph, the distance between all genes is 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 μ of the average distances of all genes 0 and standard deviation σ 0 , and then set the distance threshold to μ 0 -σ 0 ,According to the distance threshold, genes whose average distance between them and other genes is greater than the distance threshold are selected as elements of the feature selection set to form a preliminary feature selection set δ.

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

[0089] 6.1) Sort the genes in the preliminary feature set δ based on sort. After the sorting is completed, extract the feature subsets sequentially from the sorted list, and sequentially extract the first 5 elements, the first 10 elements, and the first 15 elements as the 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, and 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 mean square error MSE to determine the optimal number of clusters k 0 , MSE represents the mean square of the distance from each cell to the centroid of the cluster to which it belongs. Different k1 values ​​are selected to calculate the mean square error MSE after clustering. The k1 value with the smallest mean square error MSE after clustering is selected as the optimal cluster number k 0 , the optimal number of clusters k 0 The corresponding K-means clustering result 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 indicates the effect of cluster labels on gene expression, ε i is the error term, y i is the true cell gene expression feature vector;

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

[0099]

[0100] In the formula, 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 regression The 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] In the second aspect, a single-cell RNA sequencing data feature selection system based on deletion patterns includes:

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

[0106] Gene co-occurrence weighted graph construction module: The co-occurrence of genes in cells is calculated using the binary expression matrix based on the dropout mode and the chi-square statistic, the gene co-occurrence matrix is ​​determined, and the gene co-occurrence graph is constructed. The gene mutual exclusion matrix is ​​calculated through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout mode, the gene mutually exclusive pairs are determined, and the gene co-occurrence weighted graph is obtained;

[0107] Graph node embedding vector generation module: Based on the first-order similarity and second-order similarity 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;

[0108] Graph node importance calculation module: Based on the improved random walk algorithm, the gene co-occurrence weighted graph is walked to obtain the probability of the gene co-occurrence weighted graph node being visited, and the gene mutually exclusive pairs are used 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: Sort the preliminary feature gene selection set based on the node importance of the global information of the gene co-occurrence weighted graph, extract feature subsets of different sizes sequentially, 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.

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

[0112] Beneficial Effects

[0113] The technical effect of the present invention is that the data preprocessing and binarization module are used to remove noise and irrelevant information, thereby improving the accuracy of subsequent analysis. The binary expression matrix can highlight the dropout information, which is convenient for subsequent feature selection and co-occurrence analysis. The co-occurrence of gene pairs is calculated using the chi-square statistic, which can quantify the expression consistency of genes in different cells, providing strong evidence for revealing the cell state. The determined gene co-occurrence matrix and gene mutual exclusion matrix can reveal the interaction and potential regulatory relationship between genes, providing a new perspective for understanding cell function and disease mechanism. The LINE model is trained using the local information of the gene co-occurrence graph to generate node embedding vectors reflecting the relationship between genes, providing a basis for subsequent feature selection and cluster analysis. The defined first-order similarity and second-order similarity can capture the direct and indirect interactions between genes, and improve the accuracy and interpretability of the node embedding vector. The improved random walk algorithm can make full use of the gene mutual exclusion pair as a priori knowledge, generate the node feature importance based on global information, and provide a more comprehensive basis for feature selection. The average distance between genes in the embedded space is calculated using the Euclidean distance, which can quantify the similarity and difference between genes and provide a basis for preliminary feature selection. The established threshold can be adjusted according to 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 the clustering evaluation results, the impact of different feature combinations on cell clustering can be comprehensively evaluated to ensure the effectiveness of the selected features in revealing cell heterogeneity, disease mechanisms, and drug development, providing a scientific basis for the 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 schematic diagram based on a gene co-occurrence matrix and a mutual exclusion matrix; wherein (a) is a gene co-occurrence matrix and a gene co-occurrence map, and (b) is a gene mutual exclusion matrix and a gene mutually exclusive pair;

[0116] Figure 3The LINE model of the large network embedding model of the present invention and the first-order similarity definition and the second-order similarity definition schematic diagram; wherein the first-order proximity and the second-order proximity take the node 4 as the central node, the first-order similarity represents the similarity between the neighboring nodes 6 directly connected to the node 4, and the second-order similarity represents the similarity between the node 4 and the neighboring nodes of the directly connected node 6;

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

[0118] Figure 5 Schematic diagram of 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 schematic diagrams of comparison of 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 in the MATLAB programming language, using a publicly released human peripheral blood mononuclear cell dataset (PBMC), which is 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 selecting features 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 data of scRNA-seq sequencing, and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout mode;

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

[0123] 1.1) Data preprocessing: Gene expression matrix X of size H×U raw =X is the input, H and U are the number of rows and columns of the gene expression matrix, rows represent genes, 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 mode, N = 13714, L = 2638.

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

[0126] The 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 the 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 of 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 obtained, and the gene pairs with chi-square values ​​higher than the chi-square value threshold are determined as gene co-occurrence pairs, and the gene co-occurrence unweighted graph and edge set E are constructed;

[0131] 2.3) Calculate the weights of the node edges in the gene co-occurrence graph using 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 a 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 The transpose of , 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] In the formula, λ represents the proportionality coefficient, and its initial value is set to 1. In real data sets, 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 convert 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 in steps 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 , the matrix M N×N The gene pairs corresponding to the elements with a median value of 1 are determined as mutually exclusive gene pairs. 2.5) Using the threshold mut Binarize the mutual exclusion matrix and convert 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 usually 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 convert 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 yes, 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 , the matrix M N×N The gene pairs corresponding to the elements with a median value of 1 are determined as mutually exclusive gene pairs.

[0141] 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 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 similarity and second-order similarity 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] The 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 u 1 ,u 2 ,…,u d is the value of 1, 2, ... d dimensions corresponding to the node d-dimensional vector; d = 128;

[0147] Next, based on the defined first-order similarity, the joint probability distribution between directly connected nodes i and j in the gene co-occurrence weighted graph is constructed, 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] In the formula, 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, define the empirical probability distribution for:

[0151]

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

[0153] 3.2) The second-order similarity of the model LINE 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 shared neighbor nodes;

[0155] Secondly, initialize each node in the gene co-occurrence weighted graph as the d-dimensional embedding vector of its neighboring 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 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 defines the optimization loss function of the first-order similarity 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 defines the optimization loss function of the second-order similarity 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, and use negative samples to optimize the loss function of the second-order similarity;

[0169] This method obtains negative samples by random sampling to approximate the overall loss function, avoiding the calculation of all non-existent edges, thus significantly reducing the computational cost;

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

[0171]

[0172] In the formula, 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 = O 1 +O 2 , 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 in the next iteration based on the calculated gradient at the tth iteration:

[0175]

[0176] In the formula, α is the learning rate, u i t and u i t+1 They 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 obtained by iterative update.

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

[0178] The 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 moving 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 the 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 mutually exclusive 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 exclusive pair of the current node in the cell array exclusive_path. If the current node has no neighbor node, stop directly and go to 4.3). If there is a neighbor node, 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 mutual exclusion of the current node on the information will gradually weaken as the number of steps of the random walk increases. For each step, for each pair of nodes (i, j) in the cell array exclusive_path, the value of the decay matrix decay_matrix(i, j) is updated:

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

[0187] After updating the decay matrix decay_matrix, the probability matrix is ​​also updated:

[0188] probabilities=1-decay_matrix

[0189] 4.4) If the number of random walk rounds and the maximum number of iteration steps are not reached, return to 4.2), otherwise, the random walk stops, and 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 them. 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, and then the nodes are arranged in descending order based on the size of the node importance node_importances to obtain the sort.

[0190] E. The node embedding vector based on the local information of the gene co-occurrence weighted graph uses the Euclidean distance to calculate the average distance between genes in the embedding space, 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 vector of the local information of the gene co-occurrence weighted graph, the distance between all genes is 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 μ of the average distances of all genes 0 and standard deviation σ 0 , and then set the distance threshold to μ 0 -σ 0 ,According to the distance threshold, genes whose average distance between them and other genes is greater than the distance threshold are selected as elements of the feature selection set to form 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, extract feature subsets of different sizes in sequence, 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 the sorting is completed, extract the feature subsets sequentially from the sorted list, and sequentially extract the first 5 elements, the first 10 elements, and the first 15 elements as the 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, and 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. The goal of K-means clustering is to divide the 13714 cells into 12 clusters;

[0201] 6.3) Use mean square error MSE to determine the optimal number of clusters k 0 , MSE represents the mean square of the distance from each cell to the centroid of the cluster to which it belongs. Different k1 values ​​are selected to calculate the mean square error MSE after clustering. The k1 value with the smallest mean square error MSE after clustering is selected as the optimal cluster number k 0 , the optimal number of clusters k 0 The corresponding K-means clustering result 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 indicates the effect of cluster labels on gene expression, ε i is the error term, y i is the true cell gene expression feature vector;

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

[0209]

[0210] In the formula, 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 the other three 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 data of scRNA-seq sequencing, and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout mode;

[0217] Gene co-occurrence weighted graph construction module: The co-occurrence of genes in cells is calculated using the binary expression matrix based on the dropout mode and the chi-square statistic, the gene co-occurrence matrix is ​​determined, and the gene co-occurrence graph is constructed. The gene mutual exclusion matrix is ​​calculated through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout mode, the gene mutually exclusive pairs are determined, and the gene co-occurrence weighted graph is obtained;

[0218] Graph node embedding vector generation module: Based on the first-order similarity and second-order similarity 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;

[0219] Graph node importance calculation module: Based on the improved random walk algorithm, the gene co-occurrence weighted graph is walked to obtain the probability of the gene co-occurrence weighted graph node being visited, and the gene mutually exclusive pairs are used 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: Sort the preliminary feature gene selection set based on the node importance of the global information of the gene co-occurrence weighted graph, extract feature subsets of different sizes sequentially, 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.

[0222] It should be understood that the specific implementation process of the above-mentioned unit modules refers to the method content, and the present invention will not be specifically repeated here, and the division of the above-mentioned functional module units is only a division of logical functions. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. At the same time, the above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units.

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

[0224] For the specific implementation process, please refer to the description of the aforementioned 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 foregoing 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. Further, 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 is essentially or the part that contributes to the prior art, or all or part of the technical solution can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned readable storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, etc. Various media that can store program codes.

[0227] It should be emphasized that the examples described in the present invention are illustrative rather than restrictive, and 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 data of scRNA-seq sequencing, and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout mode; B. Calculate the co-occurrence of genes in cells using the binary expression matrix based on the dropout mode and the chi-square statistic, determine the gene co-occurrence matrix, construct a gene co-occurrence map, and calculate the gene mutual exclusion matrix through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout mode, determine the gene mutual exclusion pairs, and obtain the gene co-occurrence weighted map; C. Based on the first-order similarity and second-order similarity 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 the gene co-occurrence weighted graph node being visited, and use the gene mutually exclusive pairs as prior knowledge to obtain the importance of the graph node based on the global information of the gene co-occurrence weighted graph; E. The node embedding vector based on the local information of the gene co-occurrence weighted graph uses the Euclidean distance to calculate the average distance between genes in the embedding space, 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, extract feature subsets of different sizes in sequence, 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.

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 is the input, H and U are the number of rows and columns of the gene expression matrix, rows represent genes, columns represent cells, and genes expressed in less than 6% of cells and in more than 94% of cells are filtered out. The gene expression matrix after filtering 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 mode, respectively.

3. The method according to claim 1, characterized in that The step B comprises the following specific steps: 2.1) The chi-square statistic is used to calculate the co-occurrence between each pair of genes. Based on the binary data of the 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 weights of the node edges in the gene co-occurrence graph using 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 a 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 The transpose of , calculate the mutual exclusion matrix M N×N The average value of all elements in μ mut and standard deviation σ mut , gene mutual exclusivity threshold mut Established as: threshold mut =μ mut +λ×σ mut In the formula, λ represents the proportional coefficient, and the initial value is set to 1; 2.5) Using threshold mut Binarize the mutual exclusion matrix and convert 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, 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 ​​symmetric. If not, continue to increase the threshold. If yes, 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 as mutually exclusive gene pairs.

5. The method according to claim 3, characterized in that: The 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, the joint probability distribution between directly connected nodes i and j in the gene co-occurrence weighted graph is constructed, 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: In the formula, are the d-dimensional vector representations of node i and node j respectively, and Respectively represented as the corresponding 1, 2...d-dimensional values ​​in the d-dimensional vector of node i and node j, defining the empirical probability distribution for: Where W is the normalization factor, W = ∑ i,j∈E W(i,j); 3.2) The second-order similarity of the model LINE 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 shared neighbor nodes; Secondly, initialize each node in the gene co-occurrence weighted graph as the d-dimensional embedding vector of its neighboring 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 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 defines the optimization loss function of the first-order similarity 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 defines the optimization loss function of the second-order similarity 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, randomly select R nodes that are not directly connected to the node as negative samples, and use negative samples to optimize the loss function of the second-order similarity; The loss function of the optimized second-order similarity is: In the formula, 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: In the formula, α is the learning rate, u i t and u i t+1 They 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 obtained by iterative update.

6. The method according to claim 1, characterized in that The 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 in each round of random walk num_step, initialize an empty path array path to record the moving 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 exclusive pair of the current node in the cell array exclusive_path. If the current node has no neighbor node, stop directly and go to 4.3). If there is a neighbor node, 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 mutual exclusion of the current node on the information will gradually weaken as the number of steps of the random walk increases. For each step, for each pair of nodes (i, j) in the cell array exclusive_path, the value of the decay matrix decay_matrix(i, j) is updated: decay_matrix(i,j)=decay_matrix(i,j)×decay_factor After updating the decay matrix decay_matrix, the probability matrix is ​​also updated: probabilities=1-decay_matrix 4.4) If the number of random walk rounds and the maximum number of iteration steps are not reached, return to 4.2), otherwise, the random walk stops, and 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 them. 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, and then the nodes are arranged in descending order based on the size of the node importance node_importances to obtain the sort.

7. 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 vector of the local information of the gene co-occurrence weighted graph, the distance between all genes is 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 distance of all genes, then set the distance threshold to μ0-σ0. According to the distance threshold, select genes whose average distance between genes and other genes is greater than the distance threshold as elements of the feature selection set to form a preliminary feature selection set δ.

8. 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 the sorting is completed, extract the feature subsets sequentially from the sorted list, and sequentially extract the first 5 elements, the first 10 elements, and the first 15 elements as the feature subsets until all elements in the preliminary feature selection set are selected; 6.2) Take each feature subset as the input feature set, and 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 mean square error (MSE) after clustering. Select the k1 value with the smallest mean square error (MSE) after clustering as the optimal number of clusters k0. The optimal number of clusters k0 corresponds to the K-means clustering result as 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 square error (MSE) from the regression residuals regression : In the formula, 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.

9. 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 data of scRNA-seq sequencing, and binarize the gene expression matrix of the raw sequencing data to obtain a binary expression matrix based on the dropout mode; Gene co-occurrence weighted graph construction module: The co-occurrence of genes in cells is calculated using the binary expression matrix based on the dropout mode and the chi-square statistic, the gene co-occurrence matrix is ​​determined, and the gene co-occurrence graph is constructed. The gene mutual exclusion matrix is ​​calculated through the non-co-occurrence frequency between genes in the binary expression matrix based on the dropout mode, the gene mutually exclusive pairs are determined, and the gene co-occurrence weighted graph is obtained; Graph node embedding vector generation module: Based on the first-order similarity and second-order similarity 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; Graph node importance calculation module: Based on the improved random walk algorithm, the gene co-occurrence weighted graph is walked to obtain the probability of the gene co-occurrence weighted graph node being visited, and the gene mutually exclusive pairs are used as prior knowledge to obtain the graph node importance based on the global information of the gene co-occurrence weighted graph; 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: Sort the preliminary feature gene selection set based on the node importance of the global information of the gene co-occurrence weighted graph, extract feature subsets of different sizes sequentially, 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.

10. 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 8.

Images