Spatial transcriptome cell clustering, analysis method

Abstract

A spatial transcriptome cell clustering method, comprising the steps of: preprocessing gene expression of each cell in a spatial transcriptome; generating an adjacency matrix A according to the cell coordinates of the spatial transcriptome, obtaining a graph structure representation of the spatial transcriptome cells, representing the cell gene expression by a cell feature matrix X, inputting the adjacency matrix A and the cell feature matrix X into a trained graph convolutional neural network model DGI; the graph convolutional neural network model DGI outputs a node feature representation with spatial information; after dimension reduction and clustering algorithm processing on the node feature representation, the spatial transcriptome cell type is identified and obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of bioinformatics, and particularly relates to a spatial transcriptome cell clustering and analysis method. BACKGROUND

[0002] Spatial transcriptome technology is a major breakthrough in the field of bioinformatics in recent years. This technology can simultaneously measure the spatial position and intracellular transcriptome count of a large number of cells. The data set obtained by this method contains spatial coordinate information of a large number of cells, and each cell has high-dimensional attribute information. At present, the analysis of this kind of data is mostly based on traditional probability statistical model, which is difficult to process massive global spatial information and count data, and relies on specific model priori. Among them, cell type discovery is a necessary step for many data analysis and downstream research. SUMMARY

[0003] One of the embodiments of the present application is a spatial transcriptome cell clustering method based on a graph convolutional neural network, which preprocesses the gene expression of each cell in the spatial transcriptome.

[0004] An adjacency matrix A is generated according to the cell coordinates of the spatial transcriptome, and a graph structure representation of the spatial transcriptome cell is obtained, and the cell gene expression is represented by a cell feature matrix X.

[0005] The adjacency matrix A and the cell feature matrix X are input into a trained graph convolutional neural network model DGI.

[0006] The graph convolutional neural network model DGI outputs a node feature representation with spatial information.

[0007] After the node feature representation is processed by a dimension reduction and clustering algorithm, the cell types of the spatial transcriptome are identified and obtained.

[0008] The dimension reduction processing of the node feature representation refers to reducing the dimension of the node feature by using a PCA algorithm.

[0009] The K-means++ algorithm is used to cluster the nodes after PCA dimension reduction, so as to discover different cell types. BRIEF DESCRIPTION OF DRAWINGS

[0010] The above and other objects, features and advantages of the exemplary embodiments of the present application will be more apparent from the following detailed description read in conjunction with the accompanying drawings, in which several embodiments of the present application are shown by way of example, and in which:

[0011] Figure 1 A flowchart of a spatial transcriptome cell clustering and analysis method according to one of the embodiments of the present application.

[0012] Figure 2 As a result of an example according to an embodiment of the present application, a diagram showing GO terms associated with various types of cells and corresponding significance indicators p values is shown. DETAILED DESCRIPTION

[0013] Existing cell classification methods generally draw on the strategy of scRNA-Seq, that is, without considering the information that can be brought by spatial coordinates. Graph neural networks are an important research hotspot in the field of machine learning in recent years, and have shown strong ability to process specific structured data, and have good potential to be applied to spatial transcriptome data, thereby completing cell type discovery with spatial information.

[0014] The present application aims at the defects of existing cell clustering methods that do not fully utilize the spatial distribution information of cells, rely on specific prior knowledge, and cannot meet the problem of new cell type discovery, and proposes a cell clustering method based on graph neural network, which combines machine learning technology and knowledge in the field of spatial transcriptome, so as to use the clustering method to discover new cell subtypes.

[0015] According to one or more embodiments, a cell clustering method for spatial transcriptome data based on graph network machine learning. First, the gene expression in each cell is preprocessed to obtain node features. At the same time, the adjacency matrix is calculated according to the cell coordinates, so as to convert the spatial transcriptome data into a graph structure representation, each cell being a node on the graph, and spatially adjacent cells being connected by an edge. Further, graph convolutional neural network is applied to the graph structure data to realize feature embedding of each node. Then, dimensionality reduction and clustering are performed to finally realize cell type discovery with spatial structure. In order to verify that the method can discover cell types with different characteristics, the present application uses neighbor abundance analysis and differential expression gene analysis, Gene Otology (GO) terms analysis to test the differences in spatial distribution characteristics and biological functions of different cell clustering.

[0016] The gene expression preprocessing is obtained by the following method: first, the genes with expression are excluded according to the mean, and then the genes with low difference are excluded according to the variance, then the data collected from multiple batches is processed by Scanorama to remove batch effects, and finally the data is normalized, and the sum of the total expression of the remaining genes in each cell is 10000.

[0017] The adjacency matrix refers to: according to the spatial coordinate information of each cell, the Euclidean distance between each cell is calculated, and the adjacency matrix is constructed according to the distance, specifically including:

[0018] Step 1) For a group of cells in the same batch, the two-dimensional spatial coordinates of the i-th cell are v i =(xi , y i , the coordinates of the jth cell are v j = (x j , y j ), and the Euclidean distance between the two cells is d ij = ‖v i -v j ‖, the distance matrix of the cell group is N is the number of cells;

[0019] Step 2) Based on the distance matrix obtained above By setting a suitable threshold d thres , its adjacency matrix A0 = [a ij ] N×N , wherein:

[0020] Step 3) Further adjust the adjacency matrix, introduce a weight factor λ to realize the adjustment of the weight of the space information, and the adjusted adjacency matrix is A = λ * I + (1-λ)A0, I is a diagonal matrix: I ∈ R N×N .

[0021] The graph neural network model mainly uses the idea of Deep Graph Infomax (DGI) algorithm. The basic idea of DGI is to obtain the representation vector of the original network node through the GCN layer as an encoder. Next, the graph network structure is randomly disturbed to obtain the disturbed node representation as a negative sample. Finally, the loss function of the discriminator is optimized so that the algorithm can distinguish whether the positive and negative samples come from the original network. The goal of DGI is to maximize the mutual information between global features and local features to learn an encoder E to learn the representation of the node.

[0022] The dimension reduction method is principal component analysis (PCA) dimension reduction.

[0023] The clustering method is K-means++, which clusters all cells into K specified classes.

[0024] The neighbor abundance analysis first constructs a neighbor distribution matrix C = [c ij ] K×K , wherein c ij is the number of times that the jth cell appears in the neighbors of the ith cell. The definition of neighbor is that if there is an edge directly connecting two nodes, then the cells corresponding to the two nodes are neighbor cells. Then the neighbor distribution matrix is normalized, wherein n i is the number of cells in the ith class.

[0025] The differential expression gene analysis uses the Mann-Whitney u test to find the top 200 genes significantly highly expressed in each cluster, and the genes with a significance indicator p value less than 0.5 are removed.

[0026] The GO term analysis uses the 200 highly expressed genes screened in each cluster to upload to the GO term analysis official website to obtain a series of related biological functions from high to low.

[0027] In order to realize the above method, the related calculation and analysis system involved includes the establishment of the graph structure data of the spatial transcriptome, the establishment of the spatial transcriptome data feature embedding, dimension reduction and clustering based on DGI, and the clustering result analysis module.

[0028] The embodiment of the application proposes a new spatial transcriptome data processing method, that is,

[0029] The graph convolutional neural network technology is applied to the discovery of cell types, and an end-to-end data-driven machine learning task is completed on the clustering task, and new cell subtype discovery is achieved on multiple data sets.

[0030] According to one or more embodiments, as Figure 1 As shown in the figure, a cell clustering method for spatial transcriptome based on a graph convolutional neural network comprises the following steps:

[0031] Step 1) First, the single cell gene expression is preprocessed, and the following methods are used to obtain:

[0032] Step 1.1) Calculate the average expression of each gene in all cells, and the genes below a certain threshold will be excluded;

[0033] Step 1.2) Calculate the standard deviation of the expression of each gene in all cells, and the genes below a certain threshold will be excluded;

[0034] Step 1.3) For the spatial transcriptome data obtained from different batches, which includes cell spatial coordinates and gene expression levels in each cell, Scanorama is used to remove batch effects. Batch effects are simply experimental errors caused by the number of batches performed. These errors may originate from changes in the experimental environment, personnel, or time of the experiment. Therefore, it is necessary to process the collected data to remove batch effects. The Scanorama used here is derived from: Hie B, Bryson B, Berger B. Efficient integration of heterogeneous single-cell transcriptomes using Scanorama[J]. Nature biotechnology, 2019, 37(6):685-691.

[0035] Step 1.4) Normalize the remaining genes in each cell. Taking cell i as an example, the normalization method for its j-th gene is as follows: count ij x represents the expression levels of each gene obtained after the previous step. ij That is, the j-th dimension feature of cell i can be obtained, and the feature matrix X = [x ij ] N×N Here, 1 <= i <= N, 1 <= j <= N, and i is not equal to j. N, representing the total number of cells, is determined by the dataset. The number of cells in the dataset is N.

[0036] Step 2) Extract the coordinate information of cells in the same batch, calculate the Euclidean distance between each pair of cells, and construct an adjacency matrix based on the distance. Specifically, this includes:

[0037] Step 2.1) For cell groups within the same batch, the two-dimensional coordinates of the i-th cell are v i =(x i ,y i The coordinates of the j-th cell are v. j =(x j , y If j), then the Euclidean distance between these two cells is d. ij =‖v i -v j The distance matrix of this cell group is: L represents the number of cells;

[0038] Step 2.2) Based on the distance matrix obtained above By setting an appropriate threshold d thres Its adjacency matrix A0 = [aij ] L×L ,in:

[0039] Step 2.3) Further adjust the adjacency matrix by introducing a weight factor λ to regulate the weight of spatial information. The adjusted adjacency matrix is ​​A = λ*I + (1-λ)A0, where I is a diagonal matrix: I ∈ R N×N In this example, λ = 0.8.

[0040] Step 3) After the first two steps, the spatial transcriptome data has been transformed into a graph structure representation, where the feature matrix X represents the gene expression of the cells themselves, and the adjacency matrix A represents their spatial adjacency information. These are then input into the graph convolutional neural network model DGI to fuse the gene expression of the cells with their spatial distribution, resulting in node embeddings with spatial information. The specific method is as follows:

[0041] In the experiment, the DGI model contains four graph convolutional layers. In the l-th graph convolutional layer, node features are updated according to the following formula: Where H l and H l+1 These are the input and output of the l-th layer, W. (l) This is the weight matrix used for feature transformation in this layer. It is an adjacency matrix with self-loops added, i.e. I∈R N*N , The output of the last convolutional layer will be fed into the PReLU activation function: Where 'a' is a learnable parameter.

[0042] The goal of DGI is to learn an encoder E to learn the representation of nodes: H = E(X, A). A readout function S is also needed to extract global features s from local features H: s = S(H). For these local and global features, a discriminator D is trained in DGI to evaluate each pair (h... i The degree of matching between D(h,s). i The higher the value of s, the more global information is included in the current local features. To train this discriminator, a corruption function C is introduced to generate negative samples. In the corruption function C, the original graph structure is randomly shuffled to obtain the adjacency matrix of the negative samples. Furthermore based on And X, can obtain local features on negative samples. The final loss function is:

[0043]

[0044] By maximizing h i The mutual information between s and DGI results in the output of the nodes being the node representations after learning spatial information, while also reducing the node dimensionality.

[0045] Step 4) For the node representation of the DGI output, use PCA to further reduce the dimensionality of the node features.

[0046] Step 5) Use the K-means++ algorithm to cluster the nodes after PCA dimensionality reduction to discover different cell types. When using K-means++ clustering, the robustness of the clustering results is ensured by employing multiple random initializations of cluster centers.

[0047] Step 6) Neighbor abundance analysis: First, construct a neighbor distribution matrix C = [c ij ] K×K , where c ij Let be the number of times a cell of type j appears among the neighbors of cell of type i. A neighbor is defined as follows: if there exists an edge directly connecting two nodes, then the cells corresponding to those two nodes are neighbors. Then, the neighbor distribution matrix is ​​normalized. Where n i Let be the number of cells of type i.

[0048] Step 7) Differential gene expression analysis: This method uses the Mann-Whitney U test to find the top 200 genes that are significantly overexpressed in each cluster. Genes with a significance index p-value less than 0.5 are removed.

[0049] Step 8) GO term analysis: This method uses 200 highly expressed genes selected from each cluster and uploads them to the official GO term analysis website to obtain a series of related biological functions from high to low. Background genes are all genes in this dataset.

[0050] This embodiment conducted experiments on the MERFISH dataset, which contains three batches of cells: 645,400, and 323 cells, each providing the expression levels of 10,050 genes. Using the methods described in this embodiment, the dataset was processed and analyzed, resulting in four clusters with significant differences in biological function, which correspond well to the four different stages of the cell cycle.

[0051] First, the adjacency matrix and node features are calculated to obtain graph-structured data. This data is then input into the DGI model for training, with an epoch of 5000 and an output dimension of 256. The model's output is then used as cell features incorporating spatial information. PCA is used to further reduce the cell features to 30 dimensions, and all cells are clustered in the reduced space, with the number of clusters K set to 5.

[0052] Regarding the clustering results, since one cluster contained only two cells, this method only performed further analysis on the remaining four clusters. Spatially, all four cell clusters were mainly adjacent to cells of the same type, as shown in Table 1.

[0053] Table 1. Neighbor Abundance Matrix

[0054] Cluster label C0 C1 C3 C4 C0 0.705208 0.086957 0.151457 0.055901 C1 0.072481 0.810832 0.093190 0.023497 C3 0.127258 0.093938 0.768366 0.010438 C4 0.131019 0.066069 0.029115 0.772676

[0055] In differential gene expression analysis and GO term analysis, the four cell types exhibited different biological functions and corresponded well to the four phases of the cell cycle: G1, S, G2, and M (e.g., Figure 2 (As shown in Table 2). In contrast, existing methods cannot directly distinguish cells at different stages through clustering, and there is significant overlap in GO terms among the classification results, meaning that the clustering results obtained by previous methods do not effectively differentiate cells in a biological sense.

[0056] Table 2. Average repetition rate of the top 30 GO terms in each cluster.

[0057]

[0058] In the embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the couplings or direct couplings or communication connections shown or discussed may be indirect couplings or communication connections through some interfaces, apparatuses, or units, or they may be electrical, mechanical, or other forms of connection.

[0059] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0060] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and these modifications or substitutions should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A spatial transcriptome cell clustering method, characterized in that, Includes the following steps: (1) Preprocessing of gene expression in each cell of the spatial transcriptome, specifically including: Step 1.1: Calculate the average expression level of each gene in all cells. Genes below the first set threshold will be excluded, and genes above or equal to the first set threshold will be included in the first gene set. Step 1.2: Calculate the standard deviation of the expression level of each gene in the first gene set across all cells. Genes with an expression level below the second set threshold will be excluded, and genes with an expression level above or equal to the second set threshold will be included in the second gene set. Step 1.3: For the spatial transcriptome data obtained from different batches, Scanorama is used to remove batch effects. This spatial transcriptome data includes cell spatial coordinates and gene expression levels in each cell. Step 1.4: Based on the genes in the second gene set, normalize the genes of each cell, that is, Cell i, whose j-th gene is normalized to: ; For the expression levels of each gene, That is, the j-th dimension of the obtained cell i, from which the cell feature matrix X can be obtained. , 1<=i<=N, 1<=j<=N, i is not equal to j, N is the total number of cells; (2) Based on the cell coordinates of the spatial transcriptome, an adjacency matrix A is generated to obtain a graph representation of the spatial transcriptome cells. Cell gene expression is represented by the cell feature matrix X. The method for generating the adjacency matrix A includes: Step 2.1, for cell groups within the same batch, the two-dimensional spatial coordinates of the i-th cell are: The coordinates of the j-th cell are The Euclidean distance between these two cells is The distance matrix of this cell group is The number of cells; Step 2.2, based on the distance matrix Set threshold , thus obtaining its adjacency matrix , in, ; Step 2.3, Introduce weighting factors The adjacency matrix is ​​obtained as follows It is a diagonal matrix. ; (3) Input the adjacency matrix A and the cell feature matrix X into the trained graph convolutional neural network model DGI; The graph convolutional neural network model DGI model contains four graph convolutional layers. Within the l-th graph convolutional layer, the node features are: , in and They are the first Layer inputs and outputs, This is the weight matrix used for feature transformation in this layer. It is an adjacency matrix with self-loops added, i.e. , The input to the last convolutional layer is the PReLU activation function. , in It is a learnable parameter; (4) The graph convolutional neural network model DGI outputs node feature representations with spatial information; (5) The PCA algorithm is used to reduce the dimensionality of the node features; (6) Use the K-means++ algorithm to cluster the nodes after PCA dimensionality reduction to discover different cell types.

2. The spatial transcriptome cell clustering method according to claim 1, characterized in that, The clustering results were analyzed using neighbor abundance analysis, differentially expressed gene analysis, and Gene Ontology (GO) terms analysis to examine the differences in spatial distribution characteristics and biological functions among different cell clusters. The neighbor abundance analysis First, construct the neighbor distribution matrix. ,in Let represent the number of times a cell of type j appears among the neighbors of cell of type i. A neighbor is defined as follows: if there is an edge directly connecting two nodes, then the cells corresponding to the two nodes are neighbors of each other. Then, the neighbor distribution matrix is ​​normalized. ,in Let i be the number of cells of type i; The differentially expressed gene analysis used the Mann-Whitney U test to identify the top 200 genes that were significantly overexpressed in each cluster, and removed genes with a significance index p-value less than 0.

5. Gene Ontology (GO) terms analysis was performed by uploading 200 highly expressed genes selected from each cluster to the GOterm analysis website to obtain a series of related biological functions from high to low. Background genes were selected as all genes in this dataset.

3. A spatial transcriptome cell clustering device, characterized in that, The device includes a memory; and A processor coupled to the memory, the processor being configured to execute instructions stored in the memory to implement the method as claimed in claim 1.

4. A storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method as described in claim 1.

Images