A scRNA-seq data feature learning method based on a graph autoencoder

Abstract

The application discloses a scRNA-seq data feature learning method based on a graph autoencoder, and belongs to the technical field of biological information.The application solves the problems that a traditional autoencoder model ignores the correlation between cells in scRNA-seq data and is difficult to capture effective feature representation of the scRNA-seq data.The application provides a scRNA-seq data feature learning method which can process large-scale data sets and has high scalability from the perspective of a graph autoencoder, so as to re-learn features and reconstruct an expression matrix according to gene expression profiles of single cells, capture the correlation between cells in the scRNA-seq data and capture effective feature representation of the scRNA-seq data, and obtain a matrix which can be subjected to biological analysis such as dimension reduction and clustering, so that better clustering precision and more accurate and rich biological information can be obtained, and better downstream analysis effect can be obtained.The method can be used for obtaining scRNA-seq data features.

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 a relentless pursuit for scientists. Although genome-wide association studies have successfully identified thousands of genes associated with genetic variation, the detailed mechanisms by which these genes influence disease risk through specific cell types and cellular processes remain unknown. To fill this knowledge gap, scientists have been exploring the integration of multiple data types, particularly 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, developing effective computational methods to extract useful information from scRNA-seq data has become a hot topic 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 capture the intrinsic structure and representation of inputs, using this as a basis for analysis and research. In bioinformatics, autoencoder models are also widely used to analyze gene expression data, for example, using them to predict disease risk and discover potential disease-associated genes. However, traditional autoencoder models face the following challenges when applied to scRNA-seq data. First, the complex inter-cell relationships in scRNA-seq data are often overlooked in traditional autoencoder models. Second, the sparsity of scRNA-seq data makes it difficult for traditional autoencoder models to capture effective feature representations. Therefore, developing models tailored to the characteristics of scRNA-seq data is crucial. Summary of the Invention

[0005] The purpose of this invention is to solve the problem that traditional autoencoder models ignore the correlation between cells in scRNA-seq data and have difficulty capturing effective feature representations of scRNA-seq data, and propose 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 an scRNA-seq data matrix using the original scRNA-seq data, wherein the original scRNA-seq data is the expression data of genes in the cell;

[0008] The constructed scRNA-seq data matrix is ​​preprocessed to obtain the 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 mapping result as the input of the graph autoencoder, and output the final feature representation of the selected genes through the encoding module of the graph autoencoder;

[0012] Step 5: Perform dimensionality reduction on the final feature representation of the screened genes, and use the dimensionality reduction 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 are 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: From the constructed scRNA-seq data matrix, filter out the expression data corresponding to genes expressed in less than a0 cells. That is, for any gene, count the number of cells expressing the gene based on whether the gene is expressed in each cell. If the number of cells counted 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, that is, for any cell, count the number of genes expressed in the cell based on 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. After processing each cell separately, the preprocessed data matrix Y is obtained;

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

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

[0021] Step 1: Obtain an scRNA-seq training dataset, and process the obtained training dataset using the method of step 1 to obtain a preprocessed training data matrix;

[0022] Step 2: Obtain 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 level 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 the 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.1: 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 the K smallest elements 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 regarded as a node of the KNN graph, that is, the node set ν and edge set ε of the KNN graph are obtained, and the mapping result is obtained.

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

[0030] Step 3.1: 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 the K smallest elements 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 edges of the PKNN graph are obtained, and each gene is regarded as a node of the PKNN graph, that is, the node set ν and edge set ε of the PKNN graph are obtained, and the mapping 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 dataset whose feature importance scores are higher than the threshold;

[0035] Step (2) constructs a map based on the expression data of the genes selected in step (1) in the original scRNA-seq training dataset, and represents the feature of the i-th node in the map construction 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 features of the i-th node after transformation are 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 sequence, and the reconstructed expression data of each node is output;

[0048] Step (4), calculating the total loss based on 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), initialize the number of iterations l = 1, initialize 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 Splicing is performed, and the LeakyReLU activation function is applied to the splicing result. After normalization by the softmax function, the attention coefficient α of node i to the adjacent node j is calculated ij :

[0053]

[0054] in, Indicates transposition, ∥ indicates 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, set l=l+1 and return to step (32).

[0059] Furthermore, the total loss is calculated based on the probability and reconstructed expression data output by the graph autoencoder and 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 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 training.

[0066] Furthermore, 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.

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

[0068] From the perspective of graph autoencoders, the present invention provides a scRNAseq data feature learning method that can process large-scale data sets and is highly scalable. It relearns features and reconstructs expression matrices based on the gene expression profiles of individual cells. It can capture the correlation between cells in scRNA-seq data and capture effective feature representations of scRNA-seq data. The obtained matrix can be further subjected to biological analyses such as dimensionality reduction and clustering, which can obtain better clustering accuracy and more accurate and rich biological information, and can achieve better downstream analysis effects.

[0069] For example, experiments on 17 real datasets showed 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 during the clustering process of 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 This is a visualization result diagram of the participation in different clusters when the potential dimension of the present invention is 0 to 4 on the dataset worm_neuron_cell;

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

[0079] Figure 2b This is a visualization result of the participation in different clusters when the potential dimension of the present invention is 5 to 9 on the dataset worm_neuron_cell;

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

[0081] Figure 2d This is a visualization result of the participation in different clusters when the potential dimension of the present invention is 15 to 19 on the dataset worm_neuron_cell;

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

[0083] Figure 2f This is a visualization result of the participation in different clusters when the potential dimension of the present invention is 25-29 on the dataset worm_neuron_cell;

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

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

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

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

[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 dataset worm_neuron_cell;

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

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

[0092] Figure 3e This is the visualization result of the present invention on the dataset Pbmc_5k;

[0093] Figure 3f This 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 This 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 the measurement index;

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

[0100] Specific implementation method 1: Combination Figure 1 This embodiment describes a method for learning features from scRNA-seq data based on a graph autoencoder, the method specifically comprising the following steps:

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

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

[0103] Step 2: Obtain the features of each gene in the preprocessed data matrix Y. Use the obtained gene features and the random forest model to filter 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 obtained gene features are used as the input of the random forest model, and the random forest model outputs the feature importance score of each gene, 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 mapping result as the input of the graph autoencoder, and output the final feature representation of the selected genes through the encoding module of the graph autoencoder;

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

[0108] Density-based clustering is used to cluster features of the raw scRNA-seq data, and downstream biological analyses are performed based on the clustering results. Density clustering algorithms include, but are not limited to, the Density-Based Spatial Clustering with Applied Noise (DBSCAN) algorithm. The DBSCAN algorithm identifies core points in the data by calculating density reachability between data points. Clusters are then expanded based on these core points until all points are assigned to a cluster or marked as noise points, thereby generating clustering results for the raw data.

[0109] Specific embodiment 2: This embodiment further limits specific embodiment 1. The scRNA-seq data matrix is ​​constructed using the original scRNA-seq data, specifically:

[0110] The expression data of each gene in each cell are 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] Other steps and parameters are the same as those in the first embodiment.

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

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

[0114] Step 1: From the constructed scRNA-seq data matrix, filter out the expression data corresponding to genes expressed in less than a0 cells. That is, for any gene, count the number of cells expressing the gene based on whether the gene is expressed in each cell (if the expression level of the gene in a cell is not 0, the gene is considered to be expressed in this cell). If the number of cells counted 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, that is, for any cell, count the number of genes expressed in the cell based on 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. After processing each cell separately, the preprocessed data matrix Y is obtained;

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

[0118] Other steps and parameters are the same as those in the second embodiment.

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

[0120] Step 1: Obtain an scRNA-seq training dataset, and process the obtained training dataset using the method of step 1 to obtain a preprocessed training data matrix;

[0121] Step 2: Obtain 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 level 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 the 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 sets; where k is a positive integer; in the model, each tree is grown by randomly selecting a subset of features and a subset of samples;

[0124] The performance of the random forest model is evaluated using validation 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 impact of each gene as a feature on the performance of the entire model is calculated to assess 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., obtaining the most votes) or taking the average.

[0126] Other steps and parameters are the same as those in the third embodiment.

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

[0128] Step 3.1: Based on the original scRNA-seq data, obtain the features of each gene screened in step 2 (i.e., the vector consisting of the expression level of a gene in each cell of the original scRNA-seq data), and then calculate the Euclidean distance matrix D based on 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 the K smallest elements 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 regarded as a node of the KNN graph, that is, the node set ν and edge set ε of the KNN graph are obtained, and the mapping result is obtained.

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

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

[0133] Step 3.1: Based on the original scRNA-seq data, obtain the features of each gene screened 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 based on the obtained gene features, 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 the K smallest elements 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 edges of the PKNN graph are obtained, and each gene is regarded as a node of the PKNN graph, that is, the node set ν and edge set ε of the PKNN graph are obtained, and the mapping result is obtained.

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

[0137] Specific embodiment seven: This embodiment differs from specific embodiment four in that the training process of the graph autoencoder is as follows:

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

[0139] Step (2) constructs a map based on the expression data of the genes selected in step (1) in the original scRNA-seq training dataset, and represents the feature of the i-th node in the map construction 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 features of the i-th node after transformation are 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 sequence, and the reconstructed expression data of each node is output;

[0152] Step (4), calculating the total loss based on 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 and the trained graph autoencoder is obtained.

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

[0155] Specific embodiment eight: This embodiment differs from specific embodiment 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), initialize the number of iterations l = 1, initialize Wh is the transformed feature representation i ;

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

[0158]

[0159] in, Indicates transposition, ∥ indicates 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, set l=l+1 and return to step (32).

[0165] 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 embodiment nine: This embodiment differs from specific embodiment seven in that the total loss is calculated based on the probability and reconstructed expression data output by the graph autoencoder and 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 belonging to the edge, is the probability that the negative edge h of the inner product decoder output sample 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 training.

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

[0175] Specific embodiment ten: This embodiment differs from specific embodiment 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] 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 high-dimensional space to find an embedded representation of the data in low-dimensional space, thereby generating a reduced-dimensional data representation for subsequent clustering analysis.

[0178] In order to make the purpose, technical solutions 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 collects cells from mouse bladder tissue; the Muraro dataset is a highly comprehensive database that includes content related to diverse cell types and tissues; the pbmc_5k and pbmc_10k datasets contain data on human peripheral blood mononuclear cells, respectively; the Quake_10x_ series of datasets are single-cell sequencing data based on the 10x Genomics platform, focusing on cell types and gene expression in the bladder, limb muscles, and spleen, respectively; and the worm_neuron dataset focuses on neural cells in the nematode worm. The Quake_Smart-seq2_ series of datasets are single-cell sequencing data based on Smart-seq2 technology, focusing on cell types and gene expression in the diaphragm, limb muscles, lungs, and trachea, respectively.

[0182] Experimental environment and parameter settings

[0183] The processor is a 12th Gen Intel(R) Core(TM) i5-12500H, 2.50GHz, with 16GB of RAM and a 64-bit operating system. The software is implemented in Python using PyCharm under Windows 11. 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-Boulding index (DBI) were used to evaluate the proposed method. These indicators are widely used to evaluate model performance in unsupervised learning scenarios.

[0186] Experimental results analysis

[0187] Here, the method of the present invention was 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 datasets. Figures 2a to 2j As shown in the figure, the present invention can be used to more accurately analyze the participation of potential dimensions in different clusters, which is beneficial for subsequent downstream analysis. The visualization results of the method of the present invention on the dataset Romanov, dataset worm_neuron_cell, dataset Young, dataset Muraro, dataset Pbmc_5k, dataset Quake_10x_Spleen, dataset Adam, dataset mouse bladder, dataset Quake10x Bladder, and dataset Quake 10x Limb Muscle are shown in the figure. Figures 3a to 3j As shown in the figure, the experimental results of the silhouette coefficient (SC) and the Davidson-Boulding index (DBI) are shown in the figure. Figure 4 and Figure 5 shown.

[0192] The above examples are merely illustrative of the calculation model and process of the present invention and are not intended to limit the embodiments of the present invention. Persons skilled in the art will readily appreciate that other variations or modifications based on the above description are possible. This list of embodiments is not exhaustive; however, any obvious variations or modifications derived from the technical solution of the present invention remain within the scope of protection of the present invention.

Claims

1. A method for learning features from scRNA-seq data based on graph autoencoders, characterized in that: The method specifically comprises the following steps: Step 1: constructing an scRNA-seq data matrix using the original scRNA-seq data, wherein the original scRNA-seq data is the expression data of genes in the cell; The constructed scRNA-seq data matrix is ​​preprocessed to obtain the 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; The training process of the random forest model is: Step 1: Obtain an scRNA-seq training dataset, and process the obtained training dataset using the method of step 1 to obtain a preprocessed training data matrix; Step 2: Obtain 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 level 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 the cross-validation strategy. The random forest model outputs the feature importance score of each gene in the training dataset by voting or averaging; Step 3: Build a map based on the expression data of the genes selected in step 2 in the original scRNA-seq data; The step three is implemented using the first method or the second method: First way: Step 3.1: 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 the K smallest elements 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 edges of the KNN graph are obtained. Each gene is treated 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; Second way: Step 3.1: 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 the K smallest elements 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 edges of the PKNN graph are obtained. Each gene is treated 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; Step 4: Use the mapping result as the input of the graph autoencoder, and output the final feature representation of the selected genes through the encoding module of the graph autoencoder; Step 5: Perform dimensionality reduction on the final feature representation of the screened genes, and use the dimensionality reduction 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 scRNA-seq data matrix is ​​constructed using the original scRNA-seq data, specifically: The expression data of each gene in each cell are 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: From the constructed scRNA-seq data matrix, filter out the expression data corresponding to genes expressed in less than a0 cells. That is, for any gene, count the number of cells expressing the gene based on whether the gene is expressed in each cell. If the number of cells counted 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, that is, for any cell, count the number of genes expressed in the cell based on 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. After processing each cell separately, the preprocessed data matrix Y is obtained; Wherein, a0 and b are both natural numbers.

4. The scRNA-seq data feature learning method based on graph autoencoder according to claim 1, 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 dataset whose feature importance scores are higher than the threshold; Step (2) constructs a map based on the expression data of the genes selected in step (1) in the original scRNA-seq training dataset, and represents the feature of the i-th node in the map construction 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 features of the i-th node after transformation are 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 sequence, and the reconstructed expression data of each node is output; Step (4), calculating the total loss based on 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.

5. The scRNA-seq data feature learning method based on graph autoencoder according to claim 3, 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 Splicing is performed, and the LeakyReLU activation function is applied to the splicing result. After normalization by the softmax function, the attention coefficient α of node i to the adjacent node j is calculated ij : Among them, T represents transposition, || represents 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, set l=l+1 and return to step (32).

6. The scRNA-seq data feature learning method based on graph autoencoder according to claim 4, characterized in that The total loss is calculated based on the probability and reconstructed expression data output by the graph autoencoder and 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 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 training.

7. The method for learning scRNA-seq data features based on a 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 is performed using a dimension reduction algorithm based on manifold learning.

Images