Single-cell rna sequencing data clustering method and system based on graph autoencoder

Abstract

The application discloses a single-cell RNA sequencing data clustering method and system based on a graph autoencoder, and belongs to the fields of deep learning and biological information.The application solves the problem that existing clustering methods do not consider the relationship between cells in depth, resulting in low clustering accuracy of single-cell RNA sequencing data.The application adopts a feature autoencoder to extract data features, adopts a KNN algorithm to form a graph, and uses the graph as the input of a subsequent graph autoencoder, so that the relationship between cells is considered from the perspective of the graph instead of the traditional perspective of only considering isolated points; a GAT-based graph autoencoder is introduced to optimize the graph reconstruction effect, and the weight of a neighbor node of each graph node is adaptively calculated, so that high-order relationships between cells are abstracted.The application can be applied to single-cell RNA sequencing data clustering.

Description

Technical Field

[0001] This invention belongs to the fields of deep learning and bioinformatics, specifically relating to a method and system for clustering single-cell RNA sequencing data based on graph autoencoders. Background Technology

[0002] RNA-seq technology is an emerging high-throughput sequencing technology dedicated to data analysis. Its main task is to obtain measurement data on transcript expression levels and splicing variations in RNA samples. Single-cell RNA sequencing (scRNA-seq) is a branch of RNA-seq technology, and its main function is to test RNA data from single cells. This technology offers higher resolution and more accurate cell type identification, enabling the detection of heterogeneity and changes in cell state among different cells, thus providing a better understanding of biological processes. For example, single-cell RNA-seq technology can perform high-throughput sequencing of gene expression profiles in different cell types in the brains of Alzheimer's patients, such as neurons, oligodendrocytes, and microglia, analyzing changes in their gene expression profiles, exploring key genes in the pathogenesis of Alzheimer's disease, and providing new ideas and methods for the early diagnosis and treatment of Alzheimer's disease.

[0003] Clustering algorithms offer a convenient way to analyze single-cell RNA sequencing data using visual representations, interpreting the proportion of each cell type and subtype. Current clustering methods primarily rely on distance-based k-means clustering. However, the typical k-means algorithm may require randomization to determine centroids, potentially leading to unstable clustering results and slow convergence. To address this, an improved k-means plus algorithm was proposed, which initializes centroids through a search method, resulting in faster convergence and more stable results. However, these methods do not deeply consider the relationships between cells, failing to utilize cell graphs, particularly those assigning greater weights to closely related cells. This lack of understanding prevents the accurate revelation of heterogeneity among cell types and subtypes, resulting in low clustering accuracy for single-cell RNA sequencing data using existing methods. Summary of the Invention

[0004] The purpose of this invention is to address the problem that existing clustering methods do not take into account the relationships between cells, resulting in low clustering accuracy for single-cell RNA sequencing data. Therefore, this invention proposes a single-cell RNA sequencing data clustering method and system based on a graph autoencoder.

[0005] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0006] A clustering method for single-cell RNA sequencing data based on graph autoencoders, the method specifically includes the following steps:

[0007] Step 1: Filter and sort the single-cell RNA sequencing dataset to obtain the selected single-cell RNA sequencing data;

[0008] The screened single-cell RNA sequencing data were then normalized to obtain normalized single-cell RNA sequencing data.

[0009] Step 2: Process the normalized data based on the unsupervised feature autoencoder to obtain the feature vector output by the feature autoencoder;

[0010] Step 3: Normalize the feature vectors output by the feature autoencoder to obtain normalized feature vectors, and then convert the normalized feature vectors into the adjacency matrix of the graph.

[0011] The adjacency matrix is ​​input into a GAT-based graph autoencoder and autodecoder for reconstruction, resulting in a reconstructed adjacency matrix.

[0012] Step 4: Perform preliminary clustering on each node based on the reconstructed adjacency matrix to obtain the initial clustering results; then use the K-means clustering method to perform secondary clustering on the initial clustering results to obtain the secondary clustering results.

[0013] Clustering autoencoders were established for each type of cell obtained after secondary clustering to obtain the final cell clustering results.

[0014] Furthermore, in step one, the single-cell RNA sequencing dataset is filtered and sorted to obtain the selected single-cell RNA sequencing data; specifically:

[0015] Quality control was performed on the cells and genes in the single-cell RNA sequencing dataset to obtain the remaining cells and genes after filtering. For any remaining cell after filtering, the expression level of each remaining gene in that cell was calculated as a percentage of the total expression level of that cell. The corresponding percentages of each gene were sorted in descending order, and the top 2000 genes were selected as the RNA sequencing data of that cell after filtering. Similarly, each remaining cell after filtering was processed separately to obtain the RNA sequencing data of each remaining cell after filtering.

[0016] Furthermore, the quality control of cells is achieved based on the cell's expression characteristics, gene expression levels, and gene coverage.

[0017] Gene quality control is achieved based on gene expression levels, gene coverage, and gene detection rates.

[0018] Furthermore, the normalization operation on the screened single-cell RNA sequencing data uses a logarithmic transformation method.

[0019] Furthermore, the encoder of the feature autoencoder includes two fully connected layers, each of which performs linear transformation and nonlinear activation, and the structure of the decoder is the same as that of the encoder.

[0020] The input data x is processed by the encoder, which outputs a feature vector z. The feature vector z is then used as the input to the decoder, which reconstructs the input data x to obtain the reconstructed result.

[0021] Furthermore, the training process of the feature autoencoder is as follows:

[0022] Step 2: Process the single-cell RNA sequencing data of known categories as described in Step 1, and then divide the normalized data into several batches.

[0023] Step 22: Convert the normalized data into floating-point type, and use the converted data as input data to input the feature autoencoder in batches for encoding and decoding, to obtain the encoder output and the data reconstructed by the decoder;

[0024] Steps 2 and 3: Calculate the loss based on the reconstructed data and the input data, calculate the gradient using the backpropagation algorithm, and update the model parameters using the Adam optimizer.

[0025] Furthermore, the normalization process for the feature vector output by the feature autoencoder is performed using a binarization method.

[0026] Furthermore, the process of converting the normalized feature vectors into an adjacency matrix of a graph is achieved using the weighted KNN algorithm.

[0027] Furthermore, the adjacency matrix is ​​input into a GAT-based graph autoencoder and autodecoder for reconstruction to obtain the reconstructed adjacency matrix; specifically:

[0028] Step 3: 1. Use Dropout to regularize the adjacency matrix to obtain the processed adjacency matrix;

[0029] Step 32: Construct a graph autoencoder consisting of an input layer, a first GAT layer, a second GAT layer, a third GAT layer, and an output layer. The adjacency matrix processed in Step 31 is input into the graph autoencoder through the input layer, and then passes through the first GAT layer, the second GAT layer, and the third GAT layer in sequence.

[0030] The working process of the first GAT layer is as follows:

[0031] The attention weight 'a' between node i and node j is calculated based on the processed adjacency matrix of the input. ij :

[0032] a ij =softmax(LeakyReLU(a T [Wh i |Wh j ]))

[0033] Where LeakyReLU represents the ReLU activation function with a negative slope, Wh i |Wh j The representative will combine the two, a T This represents mapping the concatenated high-dimensional features to a real number; the softmax function is the normalization function, h j h represents the feature vector of node j. i Let W represent the feature vector of node i, and let W represent the weight matrix.

[0034] Based on attention weight a ij We perform a weighted summation of the feature vectors of node i and its neighbor node j to obtain a new feature vector for node i:

[0035]

[0036] Where, h′ i h represents the new feature vector of node i. j N represents the feature vector of node j. i Let represent the set of neighboring nodes of node i, and σ represent the activation function;

[0037] Then, the new feature vector h′ of node i i As input to the second GAT layer, the new feature vector h″ of node i is obtained after processing by the second GAT layer. i Then, the new feature vector h″ i As input to the third GAT layer, the new feature vector h″′ of node i is obtained after processing by the third GAT layer. i ;

[0038] New eigenvector h″′ i Transmitted to the graph decoder through the output layer;

[0039] Step 3: After receiving the new feature vector output by the graph autoencoder, the graph autodecoder performs Dropout regularization on the new feature vector and then calculates the inner product between nodes to obtain a new adjacency matrix.

[0040] The new adjacency matrix is ​​then normalized using an activation function to obtain the reconstructed adjacency matrix, which is the final edge weight matrix.

[0041] Furthermore, the specific process of step four is as follows:

[0042] Step 41: Treat each node in the reconstructed adjacency matrix as a community, and calculate the modularity gain ΔQ resulting from moving node i to the adjacent community j. ij :

[0043]

[0044] in, k represents the sum of the degrees of node i within the current community. i,in This represents the sum of edge weights between node i and other nodes in the current community. Let represent the degree of node i and the sum of the degrees of all nodes in the current community, m represent the sum of the weights of all edges in the network, and k represent the sum of the weights of all edges in the network. i Indicates the degree of node i;

[0045] After calculating the modularity gain brought about by moving node i to each adjacent community j, move node i to the community corresponding to the largest modularity gain.

[0046] Similarly, each node in the reconstructed adjacency matrix is ​​processed separately;

[0047] Step 42: Treat each community as a new node, recalculate the modularity gain that each node can bring by moving to an adjacent community, and then select the community with the largest modularity gain for merging.

[0048] After merging the communities, calculate the total modularity Q:

[0049]

[0050] Among them, A ij k represents the edge weight between node i and node j. j c represents the degree of node j. i c represents the community to which node i belongs. j Indicates the community to which node j belongs, when c i =c j At that time, (c i ,c j The value of ) is 1 when c i ≠c j At that time, (c i ,c j The value of ) is 0;

[0051] Step 43: Repeat the process in Step 42 until the total modularity no longer increases, then stop iterating and use the communities obtained in the last iteration as the initial clustering results.

[0052] Step 4: Process the initial clustering results using distance-based K-Means clustering to obtain the final clustering results;

[0053] Step 45: For any cell type in the clustering results of Step 44, define a clustering autoencoder based on the feature vector output by the decoder from the RNA sequencing data of each cell in that cell type. Similarly, define a clustering autoencoder for each cell type in the clustering results of Step 44.

[0054] The clustering results in step 4 are optimized using the defined full clustering autoencoder, and the final clustering results are obtained after optimization.

[0055] A single-cell RNA sequencing data clustering system based on graph autoencoder, the system comprising a data preprocessing module, a feature autoencoder module, a feature processing module, a GAT-based network module, and a clustering module;

[0056] The data preprocessing module is used to filter and sort the single-cell RNA sequencing dataset to obtain the selected single-cell RNA sequencing data; then, the selected single-cell RNA sequencing data is normalized to obtain the normalized single-cell RNA sequencing data.

[0057] The feature autoencoder module is used to process the normalized data to obtain the feature vector output by the feature autoencoder.

[0058] The feature processing module is used to normalize the feature vectors and transform the normalized feature vectors into the adjacency matrix of the graph.

[0059] A network module based on GAT is used to reconstruct the adjacency matrix to obtain the reconstructed adjacency matrix.

[0060] The clustering module is used to perform clustering based on the reconstructed adjacency matrix and obtain the clustering results.

[0061] The beneficial effects of this invention are:

[0062] This invention constructs a graph that considers the weights of inter-cell relationships and performs clustering analysis based on the cell graph. Compared with traditional distance-based clustering methods, it employs a feature autoencoder to extract data features and uses the KNN algorithm to construct a graph, which serves as the input for the subsequent graph autoencoder. This approach considers inter-cell relationships from a graph perspective, rather than simply considering isolated points. Furthermore, it introduces a GAT-based graph autoencoder to optimize graph reconstruction, adaptively calculating the weights of each graph node's neighboring nodes. This abstracts higher-order relationships between cells and applies them to cell clustering, effectively improving clustering accuracy. Attached Figure Description

[0063] Figure 1 This is a flowchart of the method of the present invention;

[0064] Figure 2a This is a heatmap of gene expression levels before feature extraction in this invention.

[0065] Figure 2b This is a heatmap of gene expression data after feature extraction according to the present invention;

[0066] Figure 3a This is a diagram showing the clustering effect achieved using the method of this invention;

[0067] Figure 3b This is a clustering result diagram of a graph autoencoder using GCN;

[0068] Figure 4a This is a schematic diagram comparing the method of the present invention with the NMI of a traditional graph autoencoder using GCN;

[0069] Figure 4b This is a schematic diagram comparing the method of the present invention with the ARI of a traditional graph autoencoder using GCN;

[0070] Figure 5 This is a schematic diagram comparing the embodiments of the present invention with the current mainstream clustering methods using ARI;

[0071] In the figure, the four bars corresponding to each dataset represent the method of this invention, MAGIC, SAVER, and DCA from left to right.

[0072] Figure 6a A diagram showing the cell population distribution in healthy individuals and Alzheimer's patients;

[0073] Figure 6b This is a comparison chart showing the results of clustering scRNA-seq data of normal individuals and Alzheimer's patients using this invention;

[0074] Figure 7 This is a comparison chart of cell types and subtype content between normal individuals and Alzheimer's patients obtained through actual analysis using this invention. Detailed Implementation

[0075] Specific Implementation Method 1: Combination Figure 1 This embodiment describes a single-cell RNA sequencing data clustering method based on a graph autoencoder, which specifically includes the following steps:

[0076] Step 1: Filter and sort the single-cell RNA sequencing dataset to obtain the selected single-cell RNA sequencing data;

[0077] The screened single-cell RNA sequencing data were then normalized to obtain normalized single-cell RNA sequencing data.

[0078] Step 2: Process the normalized data based on the unsupervised feature autoencoder to obtain the feature vector output by the feature autoencoder (i.e., the output of the decoder);

[0079] Step 3: Normalize the feature vectors output by the feature autoencoder to obtain normalized feature vectors, and then convert the normalized feature vectors into the adjacency matrix of the graph.

[0080] The adjacency matrix is ​​input into a graph autoencoder and autodecoder based on GAT for reconstruction, and the higher-order relationships between cells are abstracted to obtain the reconstructed adjacency matrix.

[0081] Step 4: Perform preliminary clustering on each node based on the reconstructed adjacency matrix to obtain the initial clustering results; then use the K-means clustering method to perform secondary clustering on the initial clustering results to obtain the secondary clustering results.

[0082] Clustering autoencoders were established for each type of cell obtained after secondary clustering to further extract features and patterns from the data and obtain the final cell clustering results.

[0083] Clustering can detect changes in cell state and obtain the proportion of each cell type and its subtype.

[0084] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that, in step one, the single-cell RNA sequencing dataset is filtered and sorted to obtain the selected single-cell RNA sequencing data; specifically:

[0085] Quality control was performed on the cells and genes in the single-cell RNA sequencing dataset to obtain the remaining cells and genes after filtering. For any remaining cell after filtering, the expression level of each remaining gene in that cell was calculated as a percentage of the total expression level of that cell. The corresponding percentages of each gene were sorted in descending order, and the top 2000 genes were selected as the RNA sequencing data of that cell after filtering. Similarly, each remaining cell after filtering was processed separately to obtain the RNA sequencing data of each remaining cell after filtering.

[0086] The other steps and parameters are the same as in Specific Implementation Method 1.

[0087] This implementation method simplifies the data, giving it a more reasonable scale and distribution.

[0088] Specific Implementation Method 3: This implementation method differs from Specific Implementation Method 1 or 2 in that the quality control of cells is based on the expression characteristics, gene expression levels, and gene coverage of the cells.

[0089] Gene quality control is achieved based on gene expression levels, gene coverage, and gene detection rates.

[0090] Cell quality control is performed using cell expression characteristics, gene expression levels, and gene coverage indicators to remove low-quality cells and ensure data accuracy and reliability. Gene quality control is performed using gene expression levels, gene coverage, and gene detection rate indicators to remove low-quality genes and ensure data accuracy and reliability.

[0091] Other steps and parameters are the same as in specific implementation method one or two.

[0092] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One to Three in that the normalization operation on the screened single-cell RNA sequencing data is performed using a logarithmic transformation.

[0093] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0094] This implementation scales the numerical values ​​to a specific range according to a certain ratio, transforming the data distribution into a form that better conforms to a normal distribution, so as to facilitate comparison and analysis between different features and improve the robustness and reliability of the model.

[0095] Specific Implementation Method 5: This implementation method differs from one of the specific implementation methods 1 to 4 in that the encoder of the feature autoencoder includes two fully connected layers, each of which performs linear transformation and nonlinear activation, and the structure of the decoder is the same as that of the encoder.

[0096] The input data x is processed by the encoder, which outputs a feature vector z. The feature vector z is then used as the input to the decoder, which reconstructs the input data x to obtain the reconstructed result.

[0097] The other steps and parameters are the same as those in one of the specific implementation methods one to four.

[0098] Specific Implementation Method Six: This implementation method differs from Specific Implementation Methods One to Five in that the training process of the feature autoencoder is as follows:

[0099] Step 2: 1. Process the single-cell RNA sequencing data of known categories as described in Step 1, and then divide the normalized data into several batches; the size of each batch can be manually adjusted.

[0100] Step 22: Take the normalized data and corresponding indexes from Step 21 from the training loader, convert the normalized data into floating-point type, and input the converted data into the feature autoencoder in batches for encoding and decoding to obtain the encoder output and the data reconstructed by the decoder.

[0101] Steps 2 and 3: Calculate the loss based on the reconstructed data and the input data, calculate the gradient using the backpropagation algorithm, and update the model parameters using the Adam optimizer.

[0102] The other steps and parameters are the same as those in one of the specific implementation methods one to five.

[0103] Specific Implementation Method Seven: This implementation method differs from Specific Implementation Methods One to Six in that the normalization processing of the feature vector output by the feature autoencoder is performed using a binarization method.

[0104] The other steps and parameters are the same as those in one of the specific implementation methods one to six.

[0105] In this implementation, normalization involves setting elements in the input data greater than the average value to 1 and elements less than or equal to the average value to 0, resulting in a binarized feature vector. This reduces data complexity, converts continuous data into discrete data, making the data easier to process and analyze, and improves the robustness of the model.

[0106] Specific Implementation Method Eight: This implementation method differs from Specific Implementation Methods One to Seven in that the weighted KNN algorithm is used to transform the normalized feature vectors into the adjacency matrix of the graph.

[0107] The other steps and parameters are the same as those in any of the specific implementation methods one to seven.

[0108] Specific Implementation Method Nine: This implementation method differs from Specific Implementation Methods One through Eight in that the adjacency matrix is ​​input into a GAT-based graph autoencoder and autodecoder for reconstruction to obtain the reconstructed adjacency matrix; specifically:

[0109] Step 3: 1. Use Dropout to regularize the adjacency matrix to obtain the processed adjacency matrix;

[0110] This is achieved by multiplying the output of each neuron by a random binary mask, the value of which is 0 or 1, with a certain probability p of being 0 and a probability of 1-p of being 1.

[0111] Step 32: Construct a graph autoencoder consisting of an input layer, a first GAT layer, a second GAT layer, a third GAT layer, and an output layer. The adjacency matrix processed in Step 31 is input into the graph autoencoder through the input layer, and then passes through the first GAT layer, the second GAT layer, and the third GAT layer in sequence.

[0112] The working process of the first GAT layer is as follows:

[0113] The attention weight 'a' between node i and node j is calculated based on the processed adjacency matrix of the input. ij :

[0114] a ij =softmax(LeakyReLU(a T [Wh i |Wh j ]))

[0115] Where LeakyReLU represents the ReLU activation function with a negative slope, Wh i |Wh j The representative will combine the two, a T This represents mapping the concatenated high-dimensional features to a real number; the softmax function is the normalization function, h j h represents the feature vector of node j. i Let W represent the feature vector of node i, and let W represent the learnable weight matrix.

[0116] Based on attention weight a ij We perform a weighted summation of the feature vectors of node i and its neighbor node j to obtain a new feature vector for node i:

[0117]

[0118] Where, h′ i h represents the new feature vector of node i. j N represents the feature vector of node j. i Let represent the set of neighboring nodes of node i, and σ represent the activation function;

[0119] Then, the new feature vector h′ of node i i As input to the second GAT layer, the new feature vector h″ of node i is obtained after processing by the second GAT layer. i Then, the new feature vector h″ i As input to the third GAT layer, the new feature vector h″′ of node i is obtained after processing by the third GAT layer. i That is, attention weights and weighted summation are performed within each GAT layer;

[0120] New eigenvector h″′i Transmitted to the graph decoder through the output layer;

[0121] Step 3: After receiving the new feature vector output by the graph autoencoder, the graph autodecoder performs Dropout regularization on the new feature vector and then calculates the inner product between nodes to obtain a new adjacency matrix.

[0122] The new adjacency matrix is ​​then normalized using an activation function to obtain the reconstructed adjacency matrix, which is the final edge weight matrix.

[0123] The other steps and parameters are the same as those in one of the specific implementation methods one to eight.

[0124] Specific Implementation Method Ten: This implementation method differs from Specific Implementation Methods One to Nine in that the specific process of step four is as follows:

[0125] Step 41: Treat each node in the reconstructed adjacency matrix as a community, and calculate the modularity gain ΔQ resulting from moving node i to the adjacent community j. ij :

[0126]

[0127] in, k represents the sum of the degrees of node i within the current community. i,in This represents the sum of edge weights between node i and other nodes in the current community. Let represent the degree of node i and the sum of the degrees of all nodes in the current community, m represent the sum of the weights of all edges in the network, and k represent the sum of the weights of all edges in the network. i Indicates the degree of node i;

[0128] After calculating the modularity gain brought about by moving node i to each adjacent community j, move node i to the community corresponding to the largest modularity gain.

[0129] Similarly, each node in the reconstructed adjacency matrix is ​​processed separately;

[0130] Step 42: Treat each community as a new node, recalculate the modularity gain that each node can bring by moving to an adjacent community, and then select the community with the largest modularity gain for merging.

[0131] After merging the communities, calculate the total modularity Q:

[0132]

[0133] Among them, A ij k represents the edge weight between node i and node j. j c represents the degree of node j. ic represents the community to which node i belongs. j Indicates the community to which node j belongs, when c i =c j At that time, (c i ,c j The value of ) is 1 when c i ≠c j At that time, (c i ,c j The value of ) is 0;

[0134] Step 43: Repeat the process in Step 42 until the total modularity no longer increases, then stop iterating and use the communities obtained in the last iteration as the initial clustering results.

[0135] Step 4: Process the initial clustering results using distance-based K-Means clustering to obtain the final clustering results;

[0136] Step 4.5 For any cell type in the clustering results of Step 4, define a clustering autoencoder based on the feature vector output by the decoder (i.e., the output in Step 2) of the RNA sequencing data of each cell in that cell type. Similarly, define a clustering autoencoder for each cell type in the clustering results of Step 4.

[0137] The clustering results in step 4 are optimized using the defined full clustering autoencoder, and the final clustering results are obtained after optimization.

[0138] Previously, Louvain and K-Means clustering were used to cluster single-cell RNA sequencing data. However, these clustering methods are based on distance or similarity and may not capture all features and patterns in the data. Therefore, defining a separate clustering autoencoder for each cluster can help further extract features and patterns from the data and more accurately assign cells to different clusters, thus better distinguishing different cell types and subtypes.

[0139] The other steps and parameters are the same as those in any of the specific implementation methods one to nine.

[0140] Specific Implementation Method 11: The single-cell RNA sequencing data clustering system based on graph autoencoder described in this implementation method includes a data preprocessing module, a feature autoencoder module, a feature processing module, a GAT-based network module, and a clustering module.

[0141] The data preprocessing module is used to filter and sort the single-cell RNA sequencing dataset to obtain the selected single-cell RNA sequencing data; then, the selected single-cell RNA sequencing data is normalized to obtain the normalized single-cell RNA sequencing data.

[0142] The feature autoencoder module is used to process the normalized data to obtain the feature vector output by the feature autoencoder.

[0143] The feature processing module is used to normalize the feature vectors and transform the normalized feature vectors into the adjacency matrix of the graph.

[0144] The network module based on GAT is used to reconstruct the adjacency matrix to obtain the reconstructed adjacency matrix; the clustering module is used to perform clustering based on the reconstructed adjacency matrix to obtain the clustering results.

[0145] The data preprocessing module works as follows:

[0146] The single-cell RNA sequencing dataset was filtered and sorted to obtain the selected single-cell RNA sequencing data:

[0147] Quality control of cells in the single-cell RNA sequencing dataset is performed based on cell expression characteristics, gene expression levels, and gene coverage. Quality control of genes in the single-cell RNA sequencing dataset is performed based on gene expression levels, gene coverage, and gene detection rate, resulting in the remaining cells and genes after filtering. For any remaining cell, the proportion of each gene's expression level to the total expression level of the remaining cell is calculated. These proportions are then sorted in descending order, and the top 2000 genes are selected as the RNA sequencing data for that cell. Similarly, each remaining cell is processed individually to obtain the RNA sequencing data for each remaining cell.

[0148] The selected single-cell RNA sequencing data were then normalized using logarithmic transformation.

[0149] The encoder of the feature autoencoder module includes two fully connected layers, each of which performs linear transformation and nonlinear activation, and the structure of the decoder is the same as that of the encoder.

[0150] The input data x is processed by the encoder, which outputs a feature vector z. The feature vector z is then used as the input to the decoder, which reconstructs the input data x to obtain the reconstructed result.

[0151] The feature processing module uses binarization to normalize the feature vectors output by the feature autoencoder, and then uses the weighted KNN algorithm to transform the normalized feature vectors into the adjacency matrix of the graph.

[0152] The working process of the network module based on GAT is as follows:

[0153] The adjacency matrix is ​​regularized using Dropout to obtain the processed adjacency matrix.

[0154] Construct a graph autoencoder consisting of an input layer, a first GAT layer, a second GAT layer, a third GAT layer, and an output layer. The adjacency matrix processed in step 3.1 is input into the graph autoencoder through the input layer, and then passes through the first GAT layer, the second GAT layer, and the third GAT layer in sequence.

[0155] The working process of the first GAT layer is as follows:

[0156] The attention weight 'a' between node i and node j is calculated based on the processed adjacency matrix of the input. ij :

[0157] a ij =softmax(LeakyReLU(a T [Wh i |Wh j ]))

[0158] Where LeakyReLU represents the ReLU activation function with a negative slope, Wh i |Wh j The representative will combine the two, a T This represents mapping the concatenated high-dimensional features to a real number; the softmax function is the normalization function, h j h represents the feature vector of node j. i Let W represent the feature vector of node i, and let W represent the weight matrix.

[0159] Based on attention weight a ij We perform a weighted summation of the feature vectors of node i and its neighbor node j to obtain a new feature vector for node i:

[0160]

[0161] Where, h′ i h represents the new feature vector of node i. j N represents the feature vector of node j. i Let represent the set of neighboring nodes of node i, and σ represent the activation function;

[0162] Then, the new feature vector h′ of node i i As input to the second GAT layer, the new feature vector h″ of node i is obtained after processing by the second GAT layer. i Then, the new feature vector h″ i As input to the third GAT layer, the new feature vector h″′ of node i is obtained after processing by the third GAT layer. i ;

[0163] New eigenvector h″′ iTransmitted to the graph decoder through the output layer;

[0164] After receiving the new feature vector output by the graph autoencoder, the graph autodecoder performs Dropout regularization on the new feature vector and then calculates the inner product between nodes to obtain a new adjacency matrix.

[0165] The new adjacency matrix is ​​then normalized using an activation function to obtain the reconstructed adjacency matrix, which is the final edge weight matrix.

[0166] The clustering module works as follows:

[0167] Step 1: Treat each node in the reconstructed adjacency matrix as a community, and calculate the modularity gain ΔQ resulting from moving node i to the adjacent community j. ij :

[0168]

[0169] in, k represents the sum of the degrees of node i within the current community. i,in This represents the sum of edge weights between node i and other nodes in the current community. Let represent the degree of node i and the sum of the degrees of all nodes in the current community, m represent the sum of the weights of all edges in the network, and k represent the sum of the weights of all edges in the network. i Indicates the degree of node i;

[0170] After calculating the modularity gain brought about by moving node i to each adjacent community j, move node i to the community corresponding to the largest modularity gain.

[0171] Similarly, each node in the reconstructed adjacency matrix is ​​processed separately;

[0172] Step 2: Treat each community as a new node, recalculate the modularity gain that each node can bring by moving to an adjacent community, and then select the community with the largest modularity gain for merging.

[0173] After merging the communities, calculate the total modularity Q:

[0174]

[0175] Among them, A ij k represents the edge weight between node i and node j. j c represents the degree of node j. i c represents the community to which node i belongs. j Indicates the community to which node j belongs, when c i =c j At that time, (c i ,c jThe value of ) is 1 when c i ≠c j At that time, (c i ,c j The value of ) is 0;

[0176] Step 3: Repeat the process of Step 2 until the total modularity no longer increases, then stop iterating and use the communities obtained in the last iteration as the initial clustering results.

[0177] Step 4: Process the initial clustering results using distance-based K-Means clustering to obtain the final clustering results;

[0178] Step 5: For any cell class in the clustering results of Step 4, define a clustering autoencoder based on the feature vector output by the decoder from the RNA sequencing data of each cell in that class. Similarly, define a clustering autoencoder for each cell class in the clustering results of Step 4.

[0179] The clustering results in step 4 are optimized using the defined full clustering autoencoder, and the final clustering results are obtained after optimization.

[0180] Example

[0181] The purpose of this invention is to provide a graph autoencoder-based clustering method for single-cell RNA sequencing data. This method can construct graphs from single-cell RNA sequencing data (scRNA-seq), including data from Alzheimer's patients, that better consider the weights of inter-cell relationships. This allows for cluster analysis based on the cell graph, exploring changes in cell types and subtypes between Alzheimer's patients and healthy individuals. Compared to traditional distance-based clustering methods, this method uses a feature autoencoder to extract data features and employs the KNN algorithm to construct a graph, which serves as input for subsequent graph autoencoders. This approach considers inter-cell relationships from a graph perspective, rather than simply considering isolated points. A GAT (Graph Attention Neural Network)-based graph autoencoder is introduced to optimize graph reconstruction, adaptively calculating the weights of neighboring nodes for each graph node. This abstracts higher-order relationships between cells and applies them to tasks such as cell clustering and cell type identification. A clustering module is employed, first performing Louvain's algorithm clustering on the reconstructed graph, and then applying a clustering autoencoder to each cluster for further analysis, resulting in better final results. Compared to traditional clustering algorithms (such as MAGIC, SAVER, and DCA), the graph autoencoder-based clustering algorithm achieved better clustering results on various single-cell RNA sequencing data. Furthermore, compared with commonly used graph autoencoders employing GCN, the GAT autoencoder used in this invention also exhibits better clustering accuracy. Finally, in practical applications, this invention was used to perform clustering analysis on single-cell RNA sequencing data from healthy individuals and patients with Alzheimer's disease, visually revealing changes in cells and cell subtypes before and after the onset of the disease, demonstrating the practical effectiveness of this invention.

[0182] This invention first preprocesses the raw scRNA-seq data by filtering out low-quality cells and genes, and also removing gene data ranked below 2000 in abundance to simplify the initial data. Subsequently, the data undergoes further feature extraction processing using a feature autoencoder. This encoder aims to learn efficient representations of cells and genes for use in subsequent autoencoders. Cell graphs are generated by learning the embeddings learned in the feature autoencoder. The cell graph is represented using a KNN graph, where nodes represent individual cells and edges represent neighbor relationships between these cells. The generated KNN graph is used as input to the graph autoencoder, which employs the GAT algorithm for information transmission, resulting in a compressed graph vector. This vector is then decoded by a graph decoder to obtain a reconstructed graph. The Louvain algorithm is used to determine the number of clusters on the reconstructed graph and assign clustering autoencoders for clustering. Figure 1 A flowchart of the clustering algorithm is provided, which includes the following steps:

[0183] Step 1: Preprocess the raw scRNA-seq data:

[0184] Table 1 shows the gene expression data formats of the training and benchmark datasets. Let M be the matrix of gene expression data. i For each line of genes, M j For each column of cells, M ij This represents the expression level of a gene in the i-th row in the j-th cell, and is usually an integer.

[0185] Table 1. Gene expression matrix format of the dataset.

[0186]

[0187] Table 2 contains basic information about the datasets, including information about one dataset of Alzheimer's patients and a normal control group, and information about four datasets used to evaluate the clustering effect.

[0188] Table 2. Detailed information about the dataset.

[0189]

[0190] 1) Alzheimer's disease sample dataset

[0191] The training data (Grubman A, Chew G, Ouyang JF, et al. As a single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[J]. Nature neuroscience, 2019, 22(12):2087-2097.) comes from a study by the Homo sapiens organization and was published on NCBI. Due to limited information on how individual cell types contribute to Alzheimer's disease (AD), the organization applied single-nuclear RNA-seq (scRNA-seq) to the entorhinal cortex of the brains of twelve individuals in an experimental and control group (six with Alzheimer's disease and six without), generating a total of 13,214 high-quality cell nuclei. The dataset represents cell-type-specific gene expression patterns, i.e., gene expression levels. Gene expression refers to the process of gene transcription into RNA, an important component of intracellular biological processes. By measuring gene expression levels, we can understand gene expression patterns in different tissues and physiological states, thereby studying gene function and regulatory mechanisms within cells. The dataset is structured as a two-dimensional table, with the horizontal axis representing gene names (e.g., MALAT1) and the vertical axis representing individual cell identifiers (e.g., AAACCTGGTAGAAAGG_AD5_AD6), identifying a single-cell sample within the dataset file. In the single-cell sequencing process, each cell is assigned a unique identifier to facilitate comparisons of differences between cells and subsequent clustering analysis. The individual cell identifier consists of two parts: AAACCTGGTAGAAAGG is the cell's barcode sequence, used to distinguish between different cells; AD5_AD6 is the sample information to which the cell belongs, used to distinguish between different samples. Analyzing the sequencing data of single cells allows us to understand gene expression differences between different cells, enabling research on cell type, cell function, and cell development.

[0192] 2) Evaluation datasets from Goolam, Klein, Semrau, and Zeisel.

[0193] The benchmark dataset used in this experiment refers to the dataset after standardization (filtering, sorting, regularization), noise removal, and outlier removal. It serves as a control group to compare and evaluate the performance of previous scRNA-seq data analysis methods. These datasets are relatively small in size after processing and come from scRNA-seq experimental datasets from different laboratories. Among them, Goolam (Goolam M, Scialdone A, Graham SJL, et al. Heterogeneity in Oct4 and Sox2 targets biases cell fate in 4-cell mouseembryos[J]. Cell, 2016, 165(1):61-74.) used single-cell RNA sequencing (scRNA-seq) technology to analyze the 4-cell stage of mouse embryos. The resulting dataset contains single-cell transcriptome data of the 4-cell stage of mice, including gene expression levels and cell type information for each cell. Klein (Klein AM, Mazutis L, Akartuna I, et al. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells[J]. Cell, 2015, 161(5): 1187-1201.) sequenced the mRNA of thousands of mouse embryonic stem and differentiated cells, obtaining sequencing results. Semrau (Semrau S, Goldmann JE, Soumillon M, et al. Dynamics of lineage commitment revealed by single-cell transcriptomics of differentiating embryonic stem cells[J]. Nature communications, 2017, 8(1): 1096.) used unbiased single-cell transcriptomics to measure the gene expression dynamics of retinoic acid-driven mESC differentiation from pluripotency to lineage commitment, obtaining scRNA-seq sequencing results. Zeisel (Zeisel A, -Manchado AB, Codeluppi S, et al. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq[J]. Science, 2015, 347(6226):1138-1142.) The dataset obtained 2000 highly expressed genes from 3005 cells using scRNA sequencing technology.

[0194] Step 1.1: The raw data is filtered and sorted to partially simplify it, resulting in a more reasonable size and distribution. A quality control filter is used to remove low-quality cells based on cell expression characteristics, gene expression levels, and gene coverage. Similarly, low-quality genes are removed based on gene expression levels, gene coverage, and gene detection rate. Finally, a sorter is used to select the top 2000 genes by their abundance.

[0195] Step 1.2: Then, normalize the data by using logarithmic scaling to scale the values ​​to a specific range, transforming the data distribution into a form that better conforms to a normal distribution. This facilitates comparison and analysis between different features, improving the robustness and reliability of the model.

[0196] Step 2: Perform feature extraction and dimensionality reduction on the preprocessed data based on an unsupervised feature autoencoder:

[0197] Step 2.1: Divide the data into several parts, each part being a batch. The size of each batch can be adjusted manually.

[0198] Step 2.2: Retrieve the data and corresponding indices from the training loader, then convert the data to floating-point type and send it to the GPU. Next, encode and decode the data to obtain the reconstructed data and the output of the encoding layer.

[0199] Step 2.3: Calculate the loss based on the reconstructed and original data, compute the gradient using the backpropagation algorithm, and update the model parameters using the Adam optimizer. At the end of each epoch, concatenate the reconstructed data and the output of the encoding layer for subsequent analysis and visualization.

[0200] Step 3: Normalize the feature vectors using binarization. Use the KNN model to transform the feature vectors into an adjacency matrix of the graph. Then, reconstruct the graph represented by the adjacency matrix using a GAT-based graph autoencoder and autodecoder to abstract higher-order relationships between cells, obtaining the reconstructed graph vectors and adjacency matrix. The specific steps are as follows:

[0201] Step 3.1: The input feature vector, after feature extraction and dimensionality reduction, is normalized using binarization. Elements greater than the mean in the input data are set to 1, and elements less than or equal to the mean are set to 0, resulting in a binarized feature vector. This process converts the data into discrete features, reducing data complexity and making continuous data easier to process and analyze, while also improving the robustness of the model.

[0202] Step 3.2: Use the KNN model to reconstruct the graph from the binarized feature vectors, obtaining the adjacency matrix of the graph as the input to the graph autoencoder.

[0203] Step 3.2: For the Edgelist graph constructed using the KNN algorithm, regularization is performed using Dropout by multiplying the output of each neuron by a random binary mask. The mask value is either 0 or 1, with a probability of 0 of probability p and a probability of 1-p of probability 1.

[0204] Step 3.4: Construct the structure of the GAT autoencoder, which mainly includes an input layer, an output layer, and three GAT layers.

[0205] Step 3.5: Calculate the neighborhood-aware softmax function, similar to the standard softmax function, but when calculating the attention weight for each node, it only considers the node's neighboring nodes, rather than all nodes in the entire graph. The attention weight 'a' for each node is then derived. ij :

[0206] a ij =softmax(LeakyReLU(a T [Wh i |Wh j ]))

[0207] Where LeakyReLU represents the ReLU activation function with a negative slope, Wh i |Wh j The representative will combine the two, a T The concatenated high-dimensional features are mapped to a real number, and the softmax function is the normalization function, h. j The feature vector of node j

[0208] Step 3.6: GAT performs a weighted summation of these weights with the feature vectors of its neighboring nodes to obtain a new feature vector for node i:

[0209]

[0210] Among them, h jLet h′ represent the feature vector of node j. i N represents the new feature vector of node i. i Let W represent the set of neighboring nodes of node i, and let W represent the learnable weight matrix. ij Let σ represent the attention weights between node i and node j, and let σ represent the activation function.

[0211] Step 3.7: In the decoder, the new feature vector output by the encoder is received and first subjected to Dropout regularization. Then, the decoder calculates the inner product between nodes to obtain an adjacency matrix. Finally, the adjacency matrix is ​​normalized through an activation function to obtain the final edge weight matrix.

[0212] The workflow of a single GAT layer can be represented by the algorithm in Table 3, and the workflow of the entire GAT autoencoder can be represented by the algorithm in Table 4:

[0213] Table 3

[0214]

[0215]

[0216] Table 4

[0217]

[0218] Step 4: Generate Louvain clustering results and the number of clusters based on the edge list of the reconstructed graph. Then, use the K-means clustering method for further clustering to distinguish different cell types and subtypes. For each cell population after two clustering operations, establish a clustering autoencoder to help further extract features and patterns from the data, and more accurately assign cells to different clusters. The specific steps are as follows:

[0219] Step 4.1: First, initialize the system, treating each node as a community. Then, perform the first iteration: For each node i, calculate the modularity gain that moving it to the adjacent community j would bring.

[0220]

[0221] in, k represents the sum of the degrees of node i within the current community. i,in This represents the sum of edge weights between node i and other nodes in the current community. This represents the sum of the degrees of node i and all nodes in the current community.

[0222] The formula for calculating modularity is:

[0223]

[0224] Among them, A ij k represents the edge weight between node i and node j. i Let m represent the degree of node i, m represent the sum of the weights of all edges in the network, and c represent the degree of node i. i Indicates the community to which node i belongs, (c i ,c j ) indicates when c i =c j The value is 1 when the two are not equal, and 0 when they are not equal.

[0225] Select the largest community j and move node i into community j.

[0226] Step 4.2: After the first iteration, proceed to the second iteration: Treat the communities from the first iteration as new nodes, recalculate the modularity gain that each node can gain by moving to an adjacent community, and select Q. ij The largest community will merge communities.

[0227] Step 4.3: Repeat the second iteration, repeatedly calculating the modularity gain, until the modularity no longer increases. This completes the generation of Louvain clustering results and the number of clusters based on the edge list.

[0228] Step 4.4: Further cluster the Louvain clustering results using distance-based K-Means clustering.

[0229] Step 4.5: While Louvain and K-Means clustering methods have been used to cluster the single-cell RNA sequencing data, these methods, based on distance or similarity, may not capture all features and patterns in the data. Therefore, defining a separate clustering autoencoder for each cluster can help further extract features and patterns from the data and more accurately assign cells to different clusters, better distinguishing different cell types and subtypes.

[0230] Figure 2a and Figure 2b The images show heatmaps of the original input data after preprocessing according to this invention, and the gene expression levels of the data after feature extraction in step 2. It can be observed that after feature extraction by the feature autoencoder, the dataset exhibits a clear trend, with the data becoming more concentrated in a specific range. The original data showed higher dispersion, with scattered feature information, which is detrimental to subsequent data analysis. However, after feature extraction by the feature autoencoder, the gene expression levels are more concentrated in certain ranges. This greater concentration demonstrates the importance of this data and its greater impact on subsequent analysis.

[0231] Figure 3a and Figure 3b The diagrams show the clustering results of the present invention and the traditional graph autoencoder using GCN. It can be seen that the clustering effect using the GCN algorithm is not as good as that of the present invention. While the cell clustering effect using the GCN algorithm can more intuitively show the distribution of each cell cluster, spatially it is roughly divided into three parts: 32cell_F, 32cell_D, and 8cell_C cluster in the upper left corner of the graph, forming a mixed cell population. Cells 2 (D), EM_4 (4), and ME_4 (4) clustered in the lower right corner, forming another cell group. Cell 32 (A) was in a separate cell group. Although the GCN algorithm divided them into three cell groups, most cells were mixed in one cell cluster, making it difficult to obtain accurate clustering results and hindering cluster analysis. However, the clustering results of the GAT algorithm used in this invention, although a small number of cells were outside the cell clusters, showed that different types of cells were distributed in different clusters, resulting in a clearer classification without the mixing of different types of cells. The spatial distribution of cells was also more objective. Furthermore, compared to the clustering results of the GCN algorithm, the GAT algorithm clustered multiple cell subpopulations, achieving better results.

[0232] This invention uses ARI and NMI as evaluation criteria for clustering results.

[0233] (1) ARI (Adjusted Rand Index)

[0234] ARI is a mainstream performance evaluation metric for clustering algorithms. It measures the similarity between the original label and the cluster label by comparing them. The ARI value ranges from -1 to 1, where 1 indicates that the clustering result is completely consistent with the true label, 0 indicates that the clustering result is consistent with the random label, and -1 indicates that the clustering result is completely opposite to the true label. The calculation of ARI first defines a confusion matrix, which is formed by combining the clustering result and the true label. The elements in the confusion matrix represent whether the data in the same cluster belong to the same category in the true label. Then, the ARI coefficient is calculated based on the confusion matrix, as shown in formula (1):

[0235]

[0236] in, This represents the number of samples belonging to the same cluster in the clustering results that also belong to the same category in the true labels, a i b represents the number of samples belonging to the i-th category in the true label. j C1 represents the number of samples belonging to the j-th cluster in the clustering results, and C2 represents the number of pairwise combinations between samples.

[0237] The advantage of the ARI metric is that it can handle any number of clusters and is insensitive to the size and shape of the clustering results. However, the ARI metric also has some disadvantages, such as being unfair in comparing different numbers of clusters and not being accurate enough in handling noisy data.

[0238] (2)NMI (Normalized Mutual Information)

[0239] Similarity Mindset (NMI) is based on the concept of mutual information in information theory. It measures the amount of information shared between two clustering results and standardizes it to eliminate the influence of different cluster sizes and labels. NMI first calculates the ratio of the number of samples in each cluster to the total number of samples, thus obtaining the probability distribution of each cluster. Then, it calculates the mutual information between two clusters, i.e., the amount of information they share, used to indicate the similarity between these clusters. Finally, NMI is calculated by standardizing the mutual information to values ​​between 0 and 1, where 1 represents a perfect match and 0 represents no match. NMI is a common clustering performance evaluation metric and is often used in conjunction with other metrics such as the silhouette coefficient to comprehensively evaluate clustering effectiveness.

[0240] NMI is calculated as shown in formula (2):

[0241]

[0242] H(T) is the entropy of the true result, calculated as shown in formula (3):

[0243]

[0244] H(C) is the entropy of the clustering result labels, calculated as shown in formula (4):

[0245]

[0246] In formula (2), MI(C,T) represents the cluster belonging to cluster C in the clustering result. i And the true labels belong to cluster category T j The number of samples, |C in formulas (3) and (4) i | represents cluster C i The number of samples in |T j | Indicates category T j The number of samples in the sample.

[0247] The use of NMI as the evaluation criterion in this invention has the following advantages:

[0248] (1) NMI is insensitive to the size of clustering results and true labels. By normalizing mutual information, it makes comparisons between different datasets and different clustering algorithms more fair and reliable.

[0249] (2) Range Consistency: The value range of NMI is [0,1], where 0 indicates that there is no correlation between the clustering result and the true label, while 1 indicates complete consistency. This makes the results of NMI easy to interpret and compare, and has intuitive meaning.

[0250] (3) Insensitivity to cluster labels: NMI does not depend on the specific value of cluster labels, but rather measures the consistency of clustering results by considering the relative distribution of clusters. Therefore, even if the cluster labels in the real labels are different, as long as the cluster structures are similar, NMI can still give a high score.

[0251] Figure 4a and Figure 4b These are schematic diagrams comparing the NMI and ARI of the method of this invention with those of traditional graph autoencoders using GCN; the horizontal axis represents the number of training layers. Figure 4a The vertical axis represents the NMI coefficient. The NMI coefficient graph shows that the clustering results using the GAT algorithm as the graph autoencoder have significantly higher NMI coefficients than those using the GCN algorithm. Furthermore, the GCN algorithm only achieves good clustering results with more than 100 iterations, while the GAT model can achieve relatively high NMI coefficients even with a lower training depth. Figure 4b For ARI, GAT also showed very good results, with an ARI coefficient of over 0.7 at fewer iteration levels, demonstrating excellent clustering performance. For GCN, it also showed good clustering performance at 40 epochs, but the average clustering value of the GAT algorithm was significantly higher than that of GCN.

[0252] This invention compares the performance of several previous cell clustering algorithms with that of the proposed algorithm. MAGIC (Model-based Analysis of Genome-wide CRISPR-Cas9 Knockout) is a CRISPR-Cas9 gene knockout-based single-cell RNA sequencing (scRNA-seq) data processing method, followed by cell clustering using clustering algorithms (such as k-means and DBSCAN). SAVER (Single-cell Analysis via Expression Recovery) is a Bayesian model-based Markov chain Monte Carlo (MCMC) method for single-cell RNA sequencing (scRNA-seq) data imputation, followed by clustering analysis using clustering algorithms (such as k-Means and DBSCAN) to divide cells into different populations. DCA (Deep Count Autoencoder) is an autoencoder-based single-cell RNA sequencing (scRNA-seq) data processing algorithm that corrects technical noise and missing values ​​in the single-cell RNA-seq data, thereby improving data quality and reliability.

[0253] Figure 5 This diagram illustrates the ARI comparison between the proposed method and the aforementioned mainstream clustering methods. It shows that for the Goolam dataset, the proposed method exhibits excellent performance with an ARI coefficient of approximately 0.6. The performance of MAGIC and DCA algorithms is similar, but not as good as the proposed method. The SAVER algorithm shows the worst performance, possibly due to the data distribution pattern and dimensionality of the Goolam dataset. For the Semrau dataset, the proposed method still performs best, MAGIC and DCA perform similarly, and SAVER performs poorly. For the Zeisel dataset, the DCA algorithm, also based on a deep learning autoencoder algorithm, performs similarly to the proposed method, while MAGIC and SAVER perform poorly. For the Klein dataset, the proposed method demonstrates better performance than other algorithms because the Klein dataset and the Goolam dataset have similar data distribution patterns; both are considered "golden datasets."

[0254] Figure 6a This is a diagram showing the cell population distribution in healthy individuals and Alzheimer's patients. Figure 6bThis image shows the distribution of cell populations after clustering in normal individuals and Alzheimer's patients using the method of this invention. This invention analyzes and clusters scRNA-seq data, resulting in 10 cell clusters, including microglia, neurons, oligodendrocyte progenitor cells (OPCs), astrocytes, and 6 oligodendrocyte subclusters (Oligo_1, Oligo_2, Oligo_3, Oligo_4, Oligo_5, Oligo_6).

[0255] Figure 7 This is a comparison chart showing the cell types and subtype content of normal individuals and Alzheimer's patients obtained after actual analysis using this invention: Among the various subtypes of oligodendrocytes, Oligo_2, Oligo_3, Oligo_4, and Oligo_5 have a higher proportion in Alzheimer's patients, while Oligo_1 and Oligo_6 have a higher proportion in the normal control group. The proportions of astrocytes and OPC cells are also higher in the normal control group than in Alzheimer's patients. It can be seen that this invention is indeed effective in revealing changes in cell types and subtypes after Alzheimer's disease lesions, and has certain practical value for exploring the pathogenesis of Alzheimer's disease and detecting early Alzheimer's lesions.

[0256] The above examples of the present invention are merely illustrative of the computational model and process of the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is impossible to exhaustively list all possible implementations here. Any obvious variations or modifications derived from the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A clustering method for single-cell RNA sequencing data based on a graph autoencoder, characterized in that, The method specifically includes the following steps: Step 1: Filter and sort the single-cell RNA sequencing dataset to obtain the selected single-cell RNA sequencing data; The screened single-cell RNA sequencing data were then normalized to obtain normalized single-cell RNA sequencing data. Step 2: Process the normalized data based on the unsupervised feature autoencoder to obtain the feature vector output by the feature autoencoder; Step 3: Normalize the feature vectors output by the feature autoencoder to obtain normalized feature vectors, and then convert the normalized feature vectors into the adjacency matrix of the graph. The adjacency matrix is ​​input into a GAT-based graph autoencoder and autodecoder for reconstruction, resulting in a reconstructed adjacency matrix. Step 4: Perform preliminary clustering on each node based on the reconstructed adjacency matrix to obtain the initial clustering results; then use the K-means clustering method to perform secondary clustering on the initial clustering results to obtain the secondary clustering results. Clustering autoencoders were established for each type of cell obtained after secondary clustering to obtain the final cell clustering results; The specific process of step four is as follows: Step 41: Treat each node in the reconstructed adjacency matrix as a community, and calculate the nodes... Move to the neighboring community Modularity gain brought about by : in, Represents a node The sum of degrees within the current community, Represents a node The sum of edge weights between the node and the node in the current community. Represents a node The degree of the node is the sum of the degrees of all nodes in the current community. This represents the sum of the weights of all edges in the network. Represents a node The degree; Calculate the nodes Moved to various neighboring communities After the modularity gain brought by the process, the nodes are then... Move to the community corresponding to the largest modularity gain; Similarly, each node in the reconstructed adjacency matrix is ​​processed separately; Step 42: Treat each community as a new node, recalculate the modularity gain that each node can bring by moving to an adjacent community, and then select the community with the largest modularity gain for merging. After merging the communities, the total modularity is then calculated. : in, Represents a node and nodes Edge weights between them Represents a node The degree, Represents a node The community to which it belongs Represents a node The community to which it belongs, when = hour, The value is 1, when ≠ hour, The value is 0; Step 43: Repeat the process in Step 42 until the total modularity no longer increases, then stop iterating and use the communities obtained in the last iteration as the initial clustering results. Step 4: Process the initial clustering results using distance-based K-Means clustering to obtain the final clustering results; Step 45: For any cell type in the clustering results of Step 44, define a clustering autoencoder based on the feature vector output by the decoder from the RNA sequencing data of each cell in that cell type. Similarly, define a clustering autoencoder for each cell type in the clustering results of Step 44. The clustering results in step 4 are optimized using the defined full clustering autoencoder, and the final clustering results are obtained after optimization.

2. The method for clustering single-cell RNA sequencing data based on graph autoencoders according to claim 1, characterized in that, In step one, the single-cell RNA sequencing dataset is filtered and sorted to obtain the selected single-cell RNA sequencing data; specifically: Quality control was performed on the cells and genes in the single-cell RNA sequencing dataset to obtain the remaining cells and genes after filtering. For any remaining cell after filtering, the expression level of each remaining gene in that cell was calculated as a percentage of the total expression level of that cell. The corresponding percentages of each gene were sorted in descending order, and the top 2000 genes were selected as the RNA sequencing data of that cell after filtering. Similarly, each remaining cell after filtering was processed separately to obtain the RNA sequencing data of each remaining cell after filtering.

3. The method for clustering single-cell RNA sequencing data based on graph autoencoders according to claim 2, characterized in that, Cell quality control is achieved based on cell expression characteristics, gene expression levels, and gene coverage. Gene quality control is achieved based on gene expression levels, gene coverage, and gene detection rates.

4. The method for clustering single-cell RNA sequencing data based on graph autoencoders according to claim 3, characterized in that, The normalization operation on the screened single-cell RNA sequencing data uses a logarithmic transformation method.

5. The method for clustering single-cell RNA sequencing data based on graph autoencoders according to claim 4, characterized in that, The encoder of the feature autoencoder includes two fully connected layers, each of which performs linear transformation and nonlinear activation, and the structure of the decoder is the same as that of the encoder. Input data After passing through the encoder, the encoder outputs the feature vector. ; eigenvectors As input to the decoder, the input data is processed by the decoder. Perform the reconstruction to obtain the reconstruction result.

6. The method for clustering single-cell RNA sequencing data based on a graph autoencoder according to claim 5, characterized in that, The normalization process for the feature vector output by the feature autoencoder is performed using a binarization method.

7. A method for clustering single-cell RNA sequencing data based on a graph autoencoder according to claim 6, characterized in that, The process of converting the normalized feature vectors into an adjacency matrix of a graph is achieved using the weighted KNN algorithm.

8. A method for clustering single-cell RNA sequencing data based on a graph autoencoder according to claim 7, characterized in that, The adjacency matrix is ​​input into a GAT-based graph autoencoder and autodecoder for reconstruction to obtain the reconstructed adjacency matrix; specifically: Step 3:

1. Use Dropout to regularize the adjacency matrix to obtain the processed adjacency matrix; Step 32: Construct a graph autoencoder consisting of an input layer, a first GAT layer, a second GAT layer, a third GAT layer, and an output layer. The adjacency matrix processed in Step 31 is input into the graph autoencoder through the input layer, and then passes through the first GAT layer, the second GAT layer, and the third GAT layer in sequence. The working process of the first GAT layer is as follows: Calculate the nodes based on the processed adjacency matrix of the input. and nodes Attention weights between : in, Indicates a negative slope Activation function The representative will combine the two. This means mapping the concatenated high-dimensional features to a real number. The function is a normalized function. Represents a node eigenvectors, Represents a node eigenvectors, Represents the weight matrix; Based on attention weight with neighboring nodes The eigenvectors are weighted and summed to obtain the node. New feature vectors: in, Represents a node The new feature vector, Represents a node eigenvectors, Represents a node The set of neighboring nodes, Indicates the activation function; Then the nodes New feature vector As input to the second GAT layer, the nodes are obtained after processing by the second GAT layer. New feature vector Then the new feature vector As input to the third GAT layer, the nodes are obtained after processing by the third GAT layer. New feature vector ; New eigenvectors Transmitted to the graph decoder through the output layer; Step 3: After receiving the new feature vector output by the graph autoencoder, the graph autodecoder performs Dropout regularization on the new feature vector and then calculates the inner product between nodes to obtain a new adjacency matrix. The new adjacency matrix is ​​then normalized using an activation function to obtain the reconstructed adjacency matrix, which is the final edge weight matrix.

9. A graph autoencoder-based single-cell RNA sequencing data clustering system, the system being used to execute the graph autoencoder-based single-cell RNA sequencing data clustering method according to any one of claims 1 to 8, characterized in that, The system includes a data preprocessing module, a feature autoencoder module, a feature processing module, a GAT-based network module, and a clustering module; The data preprocessing module is used to filter and sort the single-cell RNA sequencing dataset to obtain the selected single-cell RNA sequencing data; then, the selected single-cell RNA sequencing data is normalized to obtain the normalized single-cell RNA sequencing data. The feature autoencoder module is used to process the normalized data to obtain the feature vector output by the feature autoencoder. The feature processing module is used to normalize the feature vectors and transform the normalized feature vectors into the adjacency matrix of the graph. A network module based on GAT is used to reconstruct the adjacency matrix to obtain the reconstructed adjacency matrix. The clustering module is used to perform clustering based on the reconstructed adjacency matrix and obtain the clustering results.

Images