ScRNA-seq data feature learning method based on graph automatic encoder

Abstract

The invention discloses an scRNA-seq data feature learning method based on a graph automatic encoder, and belongs to the technical field of biological information. According to the scRNA-seq data self-encoder model, the problems that a traditional self-encoder model neglects the incidence relation between cells in the scRNA-seq data and effective feature representation of the scRNA-seq data is difficult to capture are solved. The scRNAseq data feature learning method which can process a large-scale data set and has high expandability is provided from the perspective of a graph automatic encoder, so that features are learned again according to a gene expression profile of a single cell and an expression matrix is reconstructed; according to the method, the incidence relation between cells in the scRNA-seq data and the effective feature representation of the scRNA-seq data can be captured, the obtained matrix can be subjected to biological analysis such as dimension reduction and clustering, better clustering precision and more accurate and richer biological information can be obtained, and a better downstream analysis effect can be obtained. The method disclosed by the invention can be used for acquiring scRNA-seq data characteristics.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and in particular relates to a scRNA-seq data feature learning method based on a graph autoencoder. Background Art

[0002] In the biomedical field, understanding the mechanisms of disease occurrence has always been the unremitting goal of scientists. Although genome association studies have successfully identified thousands of genes associated with genetic variation, the detailed mechanisms of how these genes affect disease risk through specific cell types and cellular processes remain unknown. In order to fill this knowledge gap, scientists have been exploring the integration of multiple data types, especially single-cell RNA sequencing data, to reveal these biological relationships.

[0003] However, due to the limitations of sequencing technology and the interference of intracellular noise, scRNA-seq data are often accompanied by significant noise and missing values, which further increases the difficulty of data processing. More importantly, most scRNA-seq data are sparse, that is, the expression levels of most genes are very low, and only a few are higher. This sparsity makes it difficult to grasp the cellular gene expression pattern as a whole. The rapid development of high-throughput technology has made single-cell RNA sequencing technology an important tool for revealing complex biological problems such as cell types, heterogeneity identification, and cell state transitions. Therefore, the development of effective computational methods to extract useful information from scRNA-seq data has become a hot spot and difficulty in current biomedical research.

[0004] In recent years, autoencoder models have achieved remarkable results in fields such as image processing. Autoencoder models can obtain the intrinsic structure and representation of the input, and use this as a basis for analysis and research. In the field of bioinformatics, autoencoder models are also widely used in the analysis of gene expression data, such as predicting disease risks and discovering potential disease-related genes through autoencoder models. However, traditional autoencoder models have the following difficulties when facing scRNA-seq data. First, there are complex correlations between cells in scRNA-seq data, which are often ignored in traditional autoencoder models. Secondly, the sparse characteristics of scRNA-seq data make it difficult for traditional autoencoder models to capture effective feature representations. Therefore, it is particularly important to develop a suitable model based on the characteristics of scRNA-seq data. Summary of the invention

[0005] The purpose of the present invention is to solve the problem that the traditional autoencoder model ignores the correlation between cells in scRNA-seq data and has difficulty in capturing the effective feature representation of scRNA-seq data, and proposes a scRNA-seq data feature learning method based on graph autoencoder.

[0006] The technical solution adopted by the present invention to solve the above technical problems is: a scRNA-seq data feature learning method based on a graph autoencoder, the method specifically comprising the following steps:

[0007] Step 1: constructing a scRNA-seq data matrix using the original scRNA-seq data, wherein the original scRNA-seq data is the expression data of genes in cells;

[0008] The constructed scRNA-seq data matrix is ​​preprocessed to obtain a preprocessed data matrix Y;

[0009] Step 2: Obtain the features of each gene in the preprocessed data matrix Y, and use the obtained gene features and the random forest model to screen out genes with feature importance scores higher than the threshold from the data matrix Y;

[0010] Step 3: Build a map based on the expression data of the genes selected in step 2 in the original scRNA-seq data;

[0011] Step 4: Use the map construction result as the input of the graph autoencoder, and output the final feature representation of the screened genes through the encoding module of the graph autoencoder;

[0012] Step 5: Perform dimensionality reduction processing on the final feature representation of the screened genes, and use the dimensionality reduction processing results as the features of the original scRNA-seq data.

[0013] Furthermore, the scRNA-seq data matrix is ​​constructed using the original scRNA-seq data, specifically:

[0014] The expression data of each gene in each cell is taken as a row of the scRNA-seq data matrix, that is, each row in the scRNA-seq data matrix corresponds to the expression data of each gene in a cell.

[0015] Furthermore, the constructed scRNA-seq data matrix is ​​preprocessed to obtain a preprocessed data matrix Y; specifically:

[0016] Step 11: Filter out the expression data corresponding to genes expressed in less than a0 cells from the constructed scRNA-seq data matrix. That is, for any gene, count the number of cells in which the gene is expressed based on whether the gene is expressed in each cell. If the number of cells obtained by counting is less than a0, delete the expression column corresponding to the gene from the constructed scRNA-seq data matrix.

[0017] Step 1 and 2: Count the number of gene expressions in each cell respectively, that is, for any cell, count the number of genes expressed in the cell according to whether each gene is expressed in the cell;

[0018] If the number of gene expressions in a cell c is less than b, the expression row corresponding to cell c is deleted from the constructed scRNA-seq data matrix, and each cell is processed separately to obtain the preprocessed data matrix Y;

[0019] Wherein, a0 and b are both natural numbers.

[0020] Furthermore, the training process of the random forest model is:

[0021] Step 1, obtaining a scRNA-seq training data set, and processing the obtained training data set using the method of step 1 to obtain a preprocessed training data matrix;

[0022] Step 2, obtaining the characteristics of each gene in the preprocessed training data matrix, that is, for any gene, the characteristic of the gene is a vector composed of the corresponding expression amount of the gene in the preprocessed training data matrix;

[0023] Step 3: Use the preprocessed training data matrix and the characteristics of each gene to train the random forest model. The training process adopts a cross-validation strategy.

[0024] The random forest model outputs the feature importance score of each gene in the training dataset by voting or averaging.

[0025] Furthermore, the specific process of step three is:

[0026] Step 3: Based on the original scRNA-seq data, obtain the characteristics of each gene screened in step 2, and then calculate the Euclidean distance matrix D based on the obtained gene characteristics, where the element in the i-th row and j-th column of the matrix D is D i,j , D i,j represents the Euclidean distance between the features of the ith gene screened and the features of the jth gene screened;

[0027] Step 32: Select K elements with the smallest values ​​from the i-th row of the matrix D, use the genes corresponding to the selected elements as the adjacent nodes of the i-th gene, and add an edge between the i-th gene and each adjacent node of the i-th gene;

[0028] After traversing each row in the matrix D, all the edges of the KNN graph are obtained, and each gene is used as a node of the KNN graph, that is, the node set ν and edge set ε of the KNN graph are obtained, and the graph construction result is obtained.

[0029] Furthermore, the specific process of step three is:

[0030] Step 3: Based on the original scRNA-seq data, obtain the characteristics of each gene screened in step 2, and then calculate the matrix P based on the obtained gene characteristics, where the element in the i-th row and j-th column of the matrix P is P i,j , P i,j represents the Pearson correlation coefficient between the features of the ith gene screened and the features of the jth gene screened;

[0031] Step 32: Select K elements with the smallest values ​​from the i-th row of the matrix P, use the genes corresponding to the selected elements as the adjacent nodes of the i-th gene, and add an edge between the i-th gene and each adjacent node of the i-th gene;

[0032] After traversing each row in the matrix P, all the edges of the PKNN graph are obtained, and each gene is used as a node of the PKNN graph, that is, the node set ν and edge set ε of the PKNN graph are obtained, and the graph construction result is obtained.

[0033] Furthermore, the training process of the graph autoencoder is:

[0034] Step (1), using the random forest model to screen out genes in the training data set whose feature importance scores are higher than the threshold;

[0035] Step (2), build a map based on the expression data of the genes selected in step (1) in the original scRNA-seq training data set, and represent the feature of the i-th node in the map building result as h i ;

[0036] Step (3), using the mapping result of step (2) as the input of the graph autoencoder, which includes an encoding module and a decoding module;

[0037] In the encoding module, the weight matrix W is first used to perform a shared linear transformation on the features of each node to obtain the transformed features of each node. The feature of the i-th node after the transformation is represented as Wh i; Then, based on the transformed feature representation of each node and the multi-head attention mechanism, the final feature representation of each node is generated;

[0038] The decoding module includes an inner product decoder and a feature reconstruction network, and the inner product decoder includes an inner product calculation unit and an activation function layer;

[0039] In the decoding module, the final feature representation of each node output by the encoding module is used as the input of the inner product decoder and the feature reconstruction network respectively;

[0040] In the inner product decoder, the inner product calculation unit is first used to calculate the inner product of any two nodes:

[0041] A ij =z i ·z j

[0042] Among them, A ij is the inner product of the final feature representation of node i and the final feature representation of node j, z i is the final feature representation of node i, z j is the final feature representation of node j;

[0043] Reusing the activation function layer and A ij Calculate the probability that there is an edge between node i and node j

[0044]

[0045] Where σ is the sigmoid function;

[0046] The feature reconstruction network includes an input layer and m fully connected layers in sequence, and the first m-1 fully connected layers are all connected to a ReLU activation function layer;

[0047] In the feature reconstruction network, the final feature representation of each node passes through each layer of the feature reconstruction network in turn, and the reconstructed expression data of each node is output;

[0048] Step (4), calculating the total loss according to the input of the graph autoencoder and the probability and reconstructed expression data output by the decoding module, and adjusting the parameters of the graph autoencoder based on the calculated total loss;

[0049] The training is stopped until the loss converges and the trained graph autoencoder is obtained.

[0050] Furthermore, the feature representation of each node after transformation and the multi-head attention mechanism are used to generate the final feature representation of each node, specifically:

[0051] Step (31), initialization iteration number l = 1, initialization Wh is the transformed feature representation i ;

[0052] Step (32), using the attention mechanism a, the features corresponding to any two nodes i and j and The splicing is performed, and the LeakyReLU activation function is applied to the splicing result. After normalization through the softmax function, the attention coefficient α of node i to the adjacent node j is calculated ij :

[0053]

[0054] in, represents transposition, ∥ represents concatenation, Represents the set of adjacent nodes of node i;

[0055] Step (33), perform weighted summation on the adjacent node features of node i, and output the feature of node i

[0056] Step (34), determine whether l=K is satisfied;

[0057] If satisfied, the feature obtained in the last iteration As the final feature representation z of node i i ;

[0058] If not satisfied, set l=l+1 and return to step (32).

[0059] Furthermore, the total loss is calculated based on the probability and reconstructed expression data of the input of the graph autoencoder and the output of the decoding module, specifically:

[0060]

[0061] in, is the total loss, is the reconstruction loss, is the mean square error loss;

[0062]

[0063] Where E is the total number of edges, is the probability that the inner product decoder outputs the true edge i as belonging to the edge, is the probability that the negative edge j of the inner product decoder output sample belongs to the edge;

[0064]

[0065] Among them, f RECON(Z) is the reconstructed expression data, and X is the input of the graph autoencoder during the training process.

[0066] Furthermore, in step 5, the dimension reduction processing of the final feature representation of the screened genes adopts a dimension reduction algorithm based on manifold learning.

[0067] The beneficial effects of the present invention are:

[0068] The present invention provides a scRNA-seq data feature learning method that can process large-scale data sets and is highly scalable from the perspective of graph autoencoders, so as to relearn the features and reconstruct the expression matrix according to the gene expression spectrum of a single cell, and can capture the correlation between cells in the scRNA-seq data and the effective feature representation of the scRNA-seq data. The obtained matrix can be further subjected to biological analysis such as dimensionality reduction and clustering, and can obtain better clustering accuracy and more accurate and rich biological information, and can obtain better downstream analysis effects.

[0069] For example, experiments on 17 real data sets show that the method of the present invention can provide more accurate scRNA-seq data clustering results and has the following advantages:

[0070] (1) The present invention solves the problem of relying on manual experience and parameter sensitivity in the screening process of highly variable genes;

[0071] (2) The present invention can more accurately capture complex patterns and information in graph structured data;

[0072] (3) In practical application scenarios, the present invention can achieve high scalability and stability;

[0073] (4) The method of the present invention can process large data sets and can address the challenges brought about by the exponential growth of data in scRNA-seq experiments;

[0074] (5) The present invention can effectively learn data representation in the process of clustering scRNA-seq data and has powerful data characterization capabilities;

[0075] (6) After clustering, there is good data visualization capability. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 It is a model framework diagram of the present invention;

[0077] Figure 2a It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 0 to 4;

[0078] In the figure, Latent dimension represents the latent dimension, and Percentage represents the percentage;

[0079] Figure 2b It is a visualization result diagram of the participation degree in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 5 to 9;

[0080] Figure 2c It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 10 to 14;

[0081] Figure 2d It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 15-19;

[0082] Figure 2e It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 20-24;

[0083] Figure 2f It is a visualization result diagram of the participation degree in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 25-29;

[0084] Figure 2g It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 30-34;

[0085] Figure 2h It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 35-39;

[0086] Figure 2i It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 40-44;

[0087] Figure 2j It is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention on the data set worm_neuron_cell is 45-49;

[0088] Figure 3a This is the visualization result of the present invention on the Romanov dataset;

[0089] Figure 3bThis is the visualization result of the present invention on the data set worm_neuron_cell;

[0090] Figure 3c This is the visualization result of the present invention on the Young dataset;

[0091] Figure 3d It is the visualization result of the present invention on the Muraro dataset;

[0092] Figure 3e It is the visualization result of the present invention on the data set Pbmc_5k;

[0093] Figure 3f It is the visualization result of the present invention on the Quake_10x_Spleen dataset;

[0094] Figure 3g This is the visualization result of the present invention on the Adam dataset;

[0095] Figure 3h It is the visualization result of the present invention on the mouse bladder dataset;

[0096] Figure 3i This is the visualization result of the present invention on the Quake 10x Bladder dataset;

[0097] Figure 3j This is the visualization result of the present invention on the Quake 10x Limb Muscle dataset;

[0098] Figure 4 It is the experimental result of the present invention using the silhouette coefficient (SC) as a measurement index;

[0099] Figure 5 This is the experimental result of the present invention using the Davidson Burghardt Index (DBI) as a measurement indicator. DETAILED DESCRIPTION

[0100] Specific implementation method 1: Combination Figure 1 The present embodiment is described as follows. The present embodiment describes a method for learning features of scRNA-seq data based on a graph autoencoder, and the method specifically comprises the following steps:

[0101] Step 1: constructing a scRNA-seq data matrix using the original scRNA-seq data, wherein the original scRNA-seq data is the expression data of genes in cells;

[0102] The constructed scRNA-seq data matrix is ​​preprocessed to obtain a preprocessed data matrix Y;

[0103] Step 2: Obtain the features of each gene in the preprocessed data matrix Y, and use the obtained gene features and the random forest model to screen out genes with feature importance scores higher than a threshold (the threshold value can be set according to actual conditions) from the data matrix Y;

[0104] The gene features to be obtained are used as the input of the random forest model, and the feature importance score of each gene is output through the random forest model, and then the genes with feature importance scores higher than the threshold are screened out;

[0105] Step 3: Build a map based on the expression data of the genes selected in step 2 in the original scRNA-seq data;

[0106] Step 4: Use the map construction result as the input of the graph autoencoder, and output the final feature representation of the screened genes through the encoding module of the graph autoencoder;

[0107] Step 5: Perform dimensionality reduction processing on the final feature representation of the screened genes, and use the dimensionality reduction processing results as the features of the original scRNA-seq data.

[0108] Based on density clustering, the features of the obtained raw scRNA-seq data are clustered, and relevant biological downstream analysis is performed based on the clustering results. Density clustering algorithms include but are not limited to density-based spatial clustering with applied noise (DBSCAN) algorithms, which identify core points in the data by calculating density accessibility between data points, and expand clustering clusters based on these core points until all points are assigned to a cluster or marked as noise points, thereby obtaining clustering results of the raw data.

[0109] Specific implementation method 2: This implementation method is a further limitation of specific implementation method 1. The method of constructing a scRNA-seq data matrix using the original scRNA-seq data is as follows:

[0110] The expression data of each gene in each cell is taken as a row of the scRNA-seq data matrix, that is, each row in the scRNA-seq data matrix corresponds to the expression data of each gene in a cell.

[0111] The other steps and parameters are the same as those in the first embodiment.

[0112] For example, the gene expression data in a cell are recorded as d1, d2, ..., d n , then all gene expression data in the cell are taken as a row of the scRNA-seq data matrix.

[0113] Specific implementation method three: This implementation method is a further limitation of specific implementation method two, wherein the constructed scRNA-seq data matrix is ​​preprocessed to obtain a preprocessed data matrix Y; specifically:

[0114] Step 1: Filter out the expression data corresponding to genes expressed in less than a0 cells from the constructed scRNA-seq data matrix. That is, for any gene, count the number of cells in which the gene is expressed based on whether the gene is expressed in each cell (if the expression of the gene in a cell is not 0, the gene is considered to be expressed in the cell). If the number of cells obtained by counting is less than a0, delete the expression column corresponding to the gene from the constructed scRNA-seq data matrix.

[0115] Step 1 and 2: Count the number of gene expressions in each cell respectively, that is, for any cell, count the number of genes expressed in the cell according to whether each gene is expressed in the cell;

[0116] If the number of gene expressions in a cell c is less than b, the expression row corresponding to cell c is deleted from the constructed scRNA-seq data matrix, and each cell is processed separately to obtain the preprocessed data matrix Y;

[0117] Wherein, a0 and b are both natural numbers.

[0118] The other steps and parameters are the same as those in the second specific implementation manner.

[0119] Specific implementation method 4: This implementation method is a further limitation of specific implementation method 3. The training process of the random forest model is as follows:

[0120] Step 1, obtaining a scRNA-seq training data set, and processing the obtained training data set using the method of step 1 to obtain a preprocessed training data matrix;

[0121] Step 2, obtaining the characteristics of each gene in the preprocessed training data matrix, that is, for any gene, the characteristic of the gene is a vector composed of the corresponding expression amount of the gene in the preprocessed training data matrix;

[0122] Step 3: Use the preprocessed training data matrix and the characteristics of each gene to train the random forest model. The training process adopts a cross-validation strategy.

[0123] Cross-validation strategy: Divide the preprocessed data into k subsets of similar size; in each cross-validation iteration, use k-1 subsets as training sets and the remaining subset as validation set; k is a positive integer; in the model, each tree is grown by randomly selecting feature subsets and sample subsets;

[0124] The performance of the random forest model is evaluated using validation set data, and the AUC (Area Under the Curve) indicator is used as the main basis for performance evaluation. At the same time, based on the Gini coefficient or other evaluation indicators, the influence of each gene as a feature on the performance of the entire model is calculated to evaluate the expression variability of each gene.

[0125] The random forest model outputs the feature importance score of each gene in the training dataset by voting (i.e., getting the most votes) or taking the average.

[0126] The other steps and parameters are the same as those in the third specific implementation manner.

[0127] Specific implementation method 5: This implementation method is a further limitation of specific implementation method 4. The specific process of step 3 is as follows:

[0128] Step 3: According to the original scRNA-seq data, the features of each gene selected in step 2 are obtained respectively (i.e., the vector composed of the expression level of a gene in each cell of the original scRNA-seq data), and then the Euclidean distance matrix D is calculated according to the obtained gene features, where the element in the i-th row and j-th column of the matrix D is D i,j , D i,j represents the Euclidean distance between the features of the ith gene screened and the features of the jth gene screened;

[0129] Step 32: Select K elements with the smallest values ​​from the i-th row of the matrix D, use the genes corresponding to the selected elements as the adjacent nodes of the i-th gene, and add an edge between the i-th gene and each adjacent node of the i-th gene;

[0130] After traversing each row in the matrix D, all the edges of the KNN graph are obtained, and each gene is used as a node of the KNN graph, that is, the node set ν and edge set ε of the KNN graph are obtained, and the graph construction result is obtained.

[0131] The other steps and parameters are the same as those in the fourth embodiment.

[0132] Specific implementation method 6: This implementation method is a further limitation of specific implementation method 4. The specific process of step 3 is as follows:

[0133] Step 3: According to the original scRNA-seq data, obtain the characteristics of each gene selected in step 2 (i.e., the vector composed of the expression level of a gene in each cell of the original scRNA-seq data), and then calculate the matrix P according to the obtained gene characteristics, where the element in the i-th row and j-th column of the matrix P is P i,j , P i,j represents the Pearson correlation coefficient between the features of the ith gene screened and the features of the jth gene screened;

[0134] Step 32: Select K elements with the smallest values ​​from the i-th row of the matrix P, use the genes corresponding to the selected elements as the adjacent nodes of the i-th gene, and add an edge between the i-th gene and each adjacent node of the i-th gene;

[0135] After traversing each row in the matrix P, all the edges of the PKNN graph are obtained, and each gene is used as a node of the PKNN graph, that is, the node set ν and edge set ε of the PKNN graph are obtained, and the graph construction result is obtained.

[0136] The other steps and parameters are the same as those in the fourth embodiment.

[0137] Specific implementation method 7: This implementation method is different from specific implementation method 4 in that the training process of the graph autoencoder is:

[0138] Step (1), using the random forest model to screen out genes in the training data set whose feature importance scores are higher than the threshold;

[0139] Step (2), build a map based on the expression data of the genes selected in step (1) in the original scRNA-seq training data set, and represent the feature of the i-th node in the map building result as h i ;

[0140] Step (3), using the mapping result of step (2) as the input of the graph autoencoder, which includes an encoding module and a decoding module;

[0141] In the encoding module, the weight matrix W is first used to perform a shared linear transformation on the features of each node to obtain the transformed features of each node. The feature of the i-th node after the transformation is represented as Wh i ; Then, based on the transformed feature representation of each node and the multi-head attention mechanism, the final feature representation of each node is generated;

[0142] The decoding module includes an inner product decoder and a feature reconstruction network, and the inner product decoder includes an inner product calculation unit and an activation function layer;

[0143] In the decoding module, the final feature representation of each node output by the encoding module is used as the input of the inner product decoder and the feature reconstruction network respectively;

[0144] In the inner product decoder, the inner product calculation unit is first used to calculate the inner product of any two nodes:

[0145] A ij =z i ·z j

[0146] Among them, A ij is the inner product of the final feature representation of node i and the final feature representation of node j, z i is the final feature representation of node i, z j is the final feature representation of node j;

[0147] Reusing the activation function layer and A ij Calculate the probability that there is an edge between node i and node j

[0148]

[0149] Where σ is the sigmoid function;

[0150] The feature reconstruction network includes an input layer and m fully connected layers in sequence, and the first m-1 fully connected layers are all connected to a ReLU activation function layer;

[0151] In the feature reconstruction network, the final feature representation of each node passes through each layer of the feature reconstruction network in turn, and the reconstructed expression data of each node is output;

[0152] Step (4), calculating the total loss according to the input of the graph autoencoder and the probability and reconstructed expression data output by the decoding module, and adjusting the parameters of the graph autoencoder based on the calculated total loss;

[0153] The training is stopped until the loss converges to obtain the trained graph autoencoder.

[0154] The other steps and parameters are the same as those in the fourth embodiment.

[0155] Specific implementation eight: This implementation differs from specific implementation seven in that the final feature representation of each node is generated based on the feature representation of each node after transformation and the multi-head attention mechanism, specifically:

[0156] Step (31), initialization iteration number l = 1, initialization Wh is the transformed feature representation i ;

[0157] Step (32), using a learnable attention mechanism a, for the features corresponding to any two nodes i and j and The splicing is performed, and the LeakyReLU activation function is applied to the splicing result. After normalization through the softmax function, the attention coefficient α of node i to the adjacent node j is calculated ij :

[0158]

[0159] in, represents transposition, ∥ represents concatenation, Represents the set of adjacent nodes of node i;

[0160] Step (33), perform weighted summation on the adjacent node features of node i, and output the feature of node i

[0161]

[0162] Step (34), determine whether l=K is satisfied;

[0163] If satisfied, the feature obtained in the last iteration As the final feature representation z of node i i ;

[0164] If not satisfied, set l=l+1 and return to step (32).

[0165] The other steps and parameters are the same as those in the seventh embodiment.

[0166] It should be noted that, for the current iteration, the data used in the current iteration process is always the node feature representation obtained in the previous iteration.

[0167] Specific implementation method 9: This implementation method is different from specific implementation method 7 in that the total loss is calculated based on the probability and reconstructed expression data of the input of the graph autoencoder and the output of the decoding module, specifically:

[0168]

[0169] in, is the total loss, is the reconstruction loss, is the mean square error loss;

[0170]

[0171] Where E is the total number of edges, is the probability that the inner product decoder outputs the true edge i as belonging to the edge, is the probability that the negative edge h sampled at the output of the inner product decoder belongs to an edge;

[0172]

[0173] Among them, f RECON (Z) is the reconstructed expression data, and X is the input of the graph autoencoder during the training process.

[0174] The other steps and parameters are the same as those in the seventh embodiment.

[0175] Specific implementation method ten: This implementation method is different from specific implementation method one in that, in step five, the dimension reduction processing for the final feature representation of the screened genes adopts a dimension reduction algorithm based on manifold learning.

[0176] The other steps and parameters are the same as those in the first embodiment.

[0177] Dimensionality reduction algorithms based on manifold learning include but are not limited to the unified manifold approximation and projection (UMAP) algorithm. The UMAP algorithm performs dimensionality reduction on data by maintaining the local manifold structure of the data in a high-dimensional space to find an embedded representation of the data in a low-dimensional space, thereby generating a reduced-dimensional data representation for subsequent clustering analysis.

[0178] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with experiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0179] Dataset Overview

[0180] To evaluate the performance of our method, we focused on relatively large datasets; we selected 17 real scRNA-seq datasets with known cell types.

[0181] The mouse_bladder dataset collects gene expression, proteomics and other bioinformatics data on mouse bladder cells; the mouse_ES dataset is for cells of mouse bladder tissue; the Muraro dataset is a highly comprehensive database that contains diverse cell types and tissue-related content; the pbmc_5k and pbmc_10k datasets contain human peripheral blood mononuclear cell data respectively; the Quake_10x_series datasets are single-cell sequencing data based on the 10x Genomics platform, focusing on the cell types and gene expression of the bladder, limb muscles and spleen respectively; the worm_neuron dataset focuses on the neural cell data of nematodes. The Quake_Smart-seq2_series datasets are single-cell sequencing data based on the Smart-seq2 technology, focusing on the cell types and gene expression of the diaphragm, limb muscles, lungs and trachea respectively.

[0182] Experimental environment and parameter settings

[0183] The processor is 12th Gen Intel(R) Core(TM) i5-12500H, 2.50GHz, the memory is 16GB RNA, and the operating system is 64-bit. The software is implemented in Python language in Pycharm environment with Windows 11 as the platform, Python version is 3.9.19, Tensorflow version is 2.16.1, and Torch version is 2.2.2. The hidden layer dimension is 256, the number of attention directions is set to 10, and the Adam optimizer is used with a learning rate of 10 -4 .

[0184] Evaluation indicators

[0185] In our experiments, the silhouette coefficient (SC), area under the curve (AUC) and Davidson-Borgin index (DBI) were used to evaluate the proposed method, which are widely used in model performance evaluation in unsupervised learning scenarios.

[0186] Experimental results analysis

[0187] Here, the method of the present invention is experimented on 17 real data sets, and the obtained relevant data are shown in Table 1.

[0188] Table 1

[0189]

[0190]

[0191] The above experimental results show that the method based on graph autoencoder (GCAN) of the present invention is a promising new method. The method of the present invention has achieved good performance on 17 real data sets. Figures 2a to 2j As shown in FIG. 1 , the present invention can be used to more accurately analyze the participation of potential dimensions in different clusters, which is beneficial to subsequent downstream analysis. The visualization results of the method of the present invention on the dataset Romanov, the dataset worm_neuron_cell, the dataset Young, the dataset Muraro, the dataset Pbmc_5k, the dataset Quake_10x_Spleen, the dataset Adam, the dataset mouse bladder, the dataset Quake10x Bladder, and the dataset Quake 10x Limb Muscle are shown in FIG. Figures 3a to 3j The experimental results of the silhouette coefficient (SC) and the Davidson-Botting index (DBI) are shown in Figure 4 and Figure 5 shown.

[0192] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. A scRNA-seq data feature learning method based on graph autoencoder, characterized in that: The method specifically comprises the following steps: Step 1: constructing a scRNA-seq data matrix using the original scRNA-seq data, wherein the original scRNA-seq data is the expression data of genes in cells; The constructed scRNA-seq data matrix is ​​preprocessed to obtain a preprocessed data matrix Y; Step 2: Obtain the features of each gene in the preprocessed data matrix Y, and use the obtained gene features and the random forest model to screen out genes with feature importance scores higher than the threshold from the data matrix Y; Step 3: Build a map based on the expression data of the genes selected in step 2 in the original scRNA-seq data; Step 4: Use the map construction result as the input of the graph autoencoder, and output the final feature representation of the screened genes through the encoding module of the graph autoencoder; Step 5: Perform dimensionality reduction processing on the final feature representation of the screened genes, and use the dimensionality reduction processing results as the features of the original scRNA-seq data.

2. The scRNA-seq data feature learning method based on graph autoencoder according to claim 1, characterized in that: The method of constructing the scRNA-seq data matrix using the original scRNA-seq data is as follows: The expression data of each gene in each cell is taken as a row of the scRNA-seq data matrix, that is, each row in the scRNA-seq data matrix corresponds to the expression data of each gene in a cell.

3. The scRNA-seq data feature learning method based on graph autoencoder according to claim 2, characterized in that: The constructed scRNA-seq data matrix is ​​preprocessed to obtain a preprocessed data matrix Y; specifically: Step 11: Filter out the expression data corresponding to genes expressed in less than a0 cells from the constructed scRNA-seq data matrix. That is, for any gene, count the number of cells in which the gene is expressed based on whether the gene is expressed in each cell. If the number of cells obtained by counting is less than a0, delete the expression column corresponding to the gene from the constructed scRNA-seq data matrix. Step 1 and 2: Count the number of gene expressions in each cell respectively, that is, for any cell, count the number of genes expressed in the cell according to whether each gene is expressed in the cell; If the number of gene expressions in a cell c is less than b, the expression row corresponding to cell c is deleted from the constructed scRNA-seq data matrix, and each cell is processed separately to obtain the preprocessed data matrix Y; Wherein, a0 and b are both natural numbers.

4. The scRNA-seq data feature learning method based on graph autoencoder according to claim 3, characterized in that: The training process of the random forest model is: Step 1, obtaining a scRNA-seq training data set, and processing the obtained training data set using the method of step 1 to obtain a preprocessed training data matrix; Step 2, obtaining the characteristics of each gene in the preprocessed training data matrix, that is, for any gene, the characteristic of the gene is a vector composed of the corresponding expression amount of the gene in the preprocessed training data matrix; Step 3: Use the preprocessed training data matrix and the characteristics of each gene to train the random forest model. The training process adopts a cross-validation strategy. The random forest model outputs the feature importance score of each gene in the training dataset by voting or averaging.

5. The scRNA-seq data feature learning method based on graph autoencoder according to claim 4, characterized in that: The specific process of step three is: Step 3: Based on the original scRNA-seq data, obtain the characteristics of each gene screened in step 2, and then calculate the Euclidean distance matrix D based on the obtained gene characteristics, where the element in the i-th row and j-th column of the matrix D is D i,j , D i,j represents the Euclidean distance between the features of the ith gene screened and the features of the jth gene screened; Step 32: Select K elements with the smallest values ​​from the i-th row of the matrix D, use the genes corresponding to the selected elements as the adjacent nodes of the i-th gene, and add an edge between the i-th gene and each adjacent node of the i-th gene; After traversing each row in the matrix D, all the edges of the KNN graph are obtained, and each gene is used as a node of the KNN graph, that is, the node set ν and edge set ε of the KNN graph are obtained, and the graph construction result is obtained.

6. The scRNA-seq data feature learning method based on graph autoencoder according to claim 4, characterized in that: The specific process of step three is: Step 3: Based on the original scRNA-seq data, obtain the characteristics of each gene screened in step 2, and then calculate the matrix P based on the obtained gene characteristics, where the element in the i-th row and j-th column of the matrix P is P i,j , P i,j represents the Pearson correlation coefficient between the features of the ith gene screened and the features of the jth gene screened; Step 32: Select K elements with the smallest values ​​from the i-th row of the matrix P, use the genes corresponding to the selected elements as the adjacent nodes of the i-th gene, and add an edge between the i-th gene and each adjacent node of the i-th gene; After traversing each row in the matrix P, all the edges of the PKNN graph are obtained, and each gene is used as a node of the PKNN graph, that is, the node set ν and edge set ε of the PKNN graph are obtained, and the graph construction result is obtained.

7. The scRNA-seq data feature learning method based on graph autoencoder according to claim 4, characterized in that: The training process of the graph autoencoder is: Step (1), using the random forest model to screen out genes in the training data set whose feature importance scores are higher than the threshold; Step (2), build a map based on the expression data of the genes selected in step (1) in the original scRNA-seq training data set, and represent the feature of the i-th node in the map building result as h i ; Step (3), using the mapping result of step (2) as the input of the graph autoencoder, which includes an encoding module and a decoding module; In the encoding module, the weight matrix W is first used to perform a shared linear transformation on the features of each node to obtain the transformed features of each node. The feature of the i-th node after the transformation is represented as Wh i ; Then, based on the transformed feature representation of each node and the multi-head attention mechanism, the final feature representation of each node is generated; The decoding module includes an inner product decoder and a feature reconstruction network, and the inner product decoder includes an inner product calculation unit and an activation function layer; In the decoding module, the final feature representation of each node output by the encoding module is used as the input of the inner product decoder and the feature reconstruction network respectively; In the inner product decoder, the inner product calculation unit is first used to calculate the inner product of any two nodes: And ij =from i ·from j Among them, A ij is the inner product of the final feature representation of node i and the final feature representation of node j, z i is the final feature representation of node i, z j is the final feature representation of node j; Reusing the activation function layer and A ij Calculate the probability that there is an edge between node i and node j Where σ is the sigmoid function; The feature reconstruction network includes an input layer and m fully connected layers in sequence, and the first m-1 fully connected layers are all connected to a ReLU activation function layer; In the feature reconstruction network, the final feature representation of each node passes through each layer of the feature reconstruction network in turn, and the reconstructed expression data of each node is output; Step (4), calculating the total loss according to the input of the graph autoencoder and the probability and reconstructed expression data output by the decoding module, and adjusting the parameters of the graph autoencoder based on the calculated total loss; The training is stopped until the loss converges and the trained graph autoencoder is obtained.

8. The scRNA-seq data feature learning method based on graph autoencoder according to claim 7, characterized in that: The final feature representation of each node is generated based on the feature representation of each node after transformation and the multi-head attention mechanism, specifically: Step (31), initialize the number of iterations l = 1, initialize Wh is the transformed feature representation i ; Step (32), using the attention mechanism a, the features corresponding to any two nodes i and j and The splicing is performed, and the LeakyReLU activation function is applied to the splicing result. After normalization through the softmax function, the attention coefficient α of node i to the adjacent node j is calculated ij : Among them, T means transposition, ∥ means splicing, Represents the set of adjacent nodes of node i; Step (33), perform weighted summation on the adjacent node features of node i, and output the feature of node i Step (34), determine whether l=K is satisfied; If satisfied, the feature obtained in the last iteration As the final feature representation z of node i i ; If not satisfied, set l=l+1 and return to step (32).

9. The scRNA-seq data feature learning method based on graph autoencoder according to claim 7, characterized in that: The total loss is calculated based on the probability and reconstructed expression data of the input of the graph autoencoder and the output of the decoding module, specifically: in, is the total loss, is the reconstruction loss, is the mean square error loss; Where E is the total number of edges, is the probability that the inner product decoder outputs the true edge i as belonging to the edge, is the probability that the negative edge j of the inner product decoder output sample belongs to the edge; Among them, f RECON (Z) is the reconstructed expression data, and X is the input of the graph autoencoder during the training process.

10. The scRNA-seq data feature learning method based on graph autoencoder according to claim 1, characterized in that: In the step 5, the dimension reduction process for the final feature representation of the screened genes adopts a dimension reduction algorithm based on manifold learning.

Images