A method for identifying cell communication based on multi-omics data

Abstract

The present invention proposes a method for identifying cell communication based on multi-omics data, which is achieved through the following steps: Step 1, obtain single-cell and spatial transcriptome datasets. Step 2, preprocess the single-cell and spatial transcriptomics datasets. Step 3, construct a deep neural network model. Step 4, use the preprocessed single-cell and spatial transcriptomics datasets to train the constructed deep neural network model. Step 5, use the trained deep neural network model to identify cell communication in the data to be tested. The present invention can achieve excellent results in identifying cell communication on multiple sets of single-cell and spatial transcriptomics datasets. By comparing with other latest cell communication identification methods, it shows better performance in multiple evaluation indicators. In addition, high stability and accuracy are achieved in the face of different levels of noise and missing spatial transcriptome datasets.

Description

Technical Field

[0001] The present invention relates to a method for cell communication, and specifically to a method for identifying cell communication based on multi-omics data, belonging to the technical field of bioinformatics. Background Art

[0002] Intercellular communication plays a crucial role in maintaining the normal functions of organisms, regulating development and differentiation, and controlling immune responses. Understanding and studying intercellular communication is of key significance for the progress of biology and medicine. Intercellular communication is usually mediated by a variety of membrane-bound factors, such as ligands, receptors, extracellular matrix (ECM), integrins, and connexins. The life of multicellular organisms depends on the coordination of cell activities, and this coordination of cell activities depends on the intercellular interactions (CCI) between different cell types and tissues in the organism. Therefore, the study of cell functions increasingly requires considering the spatial environment of each cell.

[0003] In recent years, several computational strategies based on ligand-receptor (L-R) gene pairs have been developed to identify CCI from scRNA-seq data, such as SingleCellSignalR, iTALK, CellPhoneDB, and CellChat. With the rapid development of scRNA-seq and ST-seq, the study of cell communication will be more in-depth and comprehensive. These advanced technologies allow the resolution of gene expression and spatial organization at the single-cell level, revealing the complex interactions between different cell types and subtypes. Single-cell sequencing technology can accurately identify the transcriptomic characteristics of individual cells, revealing cell heterogeneity and dynamic changes, while spatial transcriptome sequencing technology provides the spatial location information of cells in their original tissue environment, retaining the tissue structure and spatial relationships between cells. Combining these two technologies, researchers can construct a more comprehensive and accurate cell communication network to reveal how cells communicate through ligand-receptor interactions, extracellular matrix components, and other signaling pathways.

[0004] However, spatial transcriptome data often contains technical noise and errors. For example, due to variations during sample processing and sequencing, it may lead to inaccurate and missing data. Therefore, a method capable of accurately identifying intercellular communication from missing or partially missing data is needed to improve the accuracy and reliability of cell communication research. Current research methods usually only predict intercellular communication on a single tissue section, ignoring the connections between multiple tissue sections. Summary of the Invention

[0005] In order to solve the problem of identifying intercellular communication in three-dimensional space, the present invention further provides a method for identifying cell communication based on multi-omics data.

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

[0007] A method for identifying cell communication based on multi-omics data, and the method for identifying cell communication based on multi-omics data is implemented through the following steps:

[0008] S1: Obtain single-cell and spatial transcriptomics datasets, and the single-cell and spatial transcriptomics datasets include a dataset composed of scRNA-seq data and ST-seq data;

[0009] S2: Preprocess the single-cell and spatial transcriptomics datasets described in S1. For scRNA-seq data, first perform data cleaning to remove cells and genes with all-zero values as a preliminary screening step. Subsequently, standardize the screened gene expression matrix to reduce data bias and improve analysis accuracy. For ST-seq data, we obtain the corresponding spatial coordinates and then construct the corresponding cell adjacency matrix;

[0010] S3: Construct a deep neural network model, and the deep neural network model includes an encoder, a variational inference module, a decoder, and an adversarial regularization module. The encoder is responsible for encoding the input graph data into a low-dimensional latent space. The variational inference module introduces a probabilistic method for the embeddings generated by the encoder. The decoder is responsible for reconstructing the adjacency matrix of the graph or the probability of the existence of edges from the node embeddings in the latent space. The adversarial regularization module enhances the robustness of the model to input perturbations and noise;

[0011] S4: Train the deep neural network model constructed in S3 based on the single-cell and spatial transcriptomics datasets preprocessed in S2;

[0012] S5: Identify cell-to-cell communication for the test data based on the deep neural network model trained in S4.

[0013] Furthermore, the steps of training the deep neural network model described in S4 include:

[0014] S401: Alignment of multi-layer tissue sections;

[0015] S402: Construction of an adjacency matrix, and the construction of the adjacency matrix is the construction of an intercellular communication adjacency matrix or the construction of an inter-gene communication adjacency matrix;

[0016] S403: Construction of a cell communication model, combine GCN and VAE to construct a cell communication model VGAE-CCI based on VGAE.

[0017] Furthermore, the alignment method described in S401 is the PASTE algorithm, and the alignment method steps of the PASTE algorithm include:

[0018] S401.1: Extract feature vectors from each slice, constructed using the gene expression data or coordinates of each cell;

[0019] S401.2: Calculate the distance between each pair of cells;

[0020] S401.3: Align the slices by optimizing the alignment objective function, and calculate the distance matrix between the points for each slice;

[0021] S401.4: Optimize the rotation and scaling parameters to make the alignment result more accurate.

[0022] Furthermore, the calculation method of the intercellular distance described in S401.2 is the Euclidean distance. Assume that xi and xj are the feature vectors of two cells, then the formula for calculating their Euclidean distance dij is as follows:

[0023]

[0024] where, x ik and x jk are respectively the k-th components in the feature vector.

[0025] Furthermore, each slice calculation step in S401.3 includes:

[0026] a. Initialize the matching matrix P, defined as follows:

[0027]

[0028] where, n and m are respectively the number of points in the two slices;

[0029] b. For each slice, calculate the distance matrix between its points, and the expression is as follows:

[0030] D X (i,j) = ||X i - X j ||

[0031] D Y (k,l) = ||Y k - Y l ||

[0032] where, D X (i,j) and D Y (k,l) are respectively the distance matrices of the two cell graphs, and X and Y are respectively the point sets of the two slices;

[0033] c. Use the Sinkhorn-Knopp algorithm for iterative optimization to update the matching matrix P, and the expression is as follows:

[0034] P ← diag(u)·P·diag(v)

[0035] Wherein, u and v are vectors for normalization, and P is updated iteratively to satisfy the row and column sum constraint conditions.

[0036] Furthermore, the optimized alignment objective function expression in S401.3 is as follows:

[0037] C(i,k,j,l) = (D X (i,j) - D y (k,l)) 2

[0038]

[0039] Wherein, the GW distance is a measure of the difference between the distance matrices between cells, thereby minimizing the GW distance between two cell graphs.

[0040] Furthermore, the method for optimizing the rotation and scaling parameters in S401.4 is the generalized Procrustes analysis, and its expression is as follows:

[0041]

[0042] Wherein, R is the rotation matrix and t is the translation vector.

[0043] Furthermore, the construction steps of the cell-level communication adjacency matrix described in S402 are as follows:

[0044] a. Calculate the distance matrix according to the coordinate data of the cells. The calculation method is the Euclidean distance, and its calculation formula is as follows:

[0045]

[0046] Wherein, D (i,j) is the distance between cell i and cell j, and x ik and x jk are the coordinates of cell i and cell j in the k-th dimension;

[0047] b. Determine the nearest neighbor cells of each cell. For each cell i, find its n nearest neighbor cells and store these neighbor cells in N i . The expression of N i is as follows:

[0048] N i = {j|D ij is the first n smallest elements in D i}

[0049] Wherein, D i represents the i-th row of the distance matrix D

[0050] c. Construct an adjacency matrix, the expression is as follows:

[0051]

[0052] d. Symmetrize the constructed adjacency matrix to make it a symmetric matrix.

[0053] Furthermore, the construction steps of the gene-level communication adjacency matrix described in S402 are as follows:

[0054] a. Given gene expression data E = {e1, e2,..., eu}, where ei is the identifier of the i-th gene and u is the total number of genes,

[0055] b. Construct an adjacency matrix, the expression is as follows:

[0056]

[0057] Among them, if genes i and j in the given gene expression matrix exist in the database and correspond to each other, then ei and ej are considered related;

[0058] c. Symmetrize the constructed adjacency matrix to make it a symmetric matrix.

[0059] Furthermore, the construction steps of the cell communication model described in S403 are as follows:

[0060] a. Combine GCN and VAE to construct an end-to-end model to simultaneously learn graph structure information and latent space representation from the cell communication network. From the output of the encoder, the model learns to generate the distribution of the latent variable z, and the expression is as follows:

[0061]

[0062] Among them, x is the input data, represents a normal distribution, and are the mean and standard deviation learned through the neural network parameters and represents the posterior distribution that describes the distribution of the latent variable z given the input x;

[0063] b. Introduce a probabilistic method to use the reparameterization trick to learn the mean and variance of the latent distribution so that the gradient can be passed through the latent variable z, and the expression is as follows:

[0064] z = μ φ (x) + σ φ (x) ⊙ ∈

[0065] Among them, is the noise sampled from the standard normal distribution;

[0066] c. Use the inner product decoder to calculate the reconstructed adjacency matrix, and the expression is as follows:

[0067]

[0068] where, is the sum of the identity matrix and the relationship matrix of each node, σ is the sigmoid function, Z is the latent space representation, and Z is constrained by the adversarial regularization module from the prior Gaussian distribution. Z T is the transpose of Z.

[0069] Furthermore, the loss function of the cell communication model VGAE-CCI is composed of the reconstruction loss and the KL divergence loss combined, and is defined as follows:

[0070]

[0071] where, β is a weight parameter used to balance the reconstruction loss and the KL divergence loss, and the default setting is 1.

[0072] Furthermore, the reconstruction loss is used to measure the difference between the generated data and the original data. The reconstruction loss is the cross-entropy loss, and the expression is as follows:

[0073]

[0074] where, logp θ (x|z) is the logarithm of the likelihood function, indicating the probability of reconstructing the input x given the latent variable z.

[0075] Furthermore, the KL divergence loss is used to measure the difference between the approximate posterior distribution and the prior distribution p(z), and the expression is as follows:

[0076]

[0077] Furthermore, the propagation method between layers in the GCN neural network layer is:

[0078]

[0079] where, is the degree matrix of, and the formula is H is the node feature matrix of each layer, and W is the learnable weight matrix of each layer.

[0080] Furthermore, the potential space representation method extracts features through a convolutional layer and generates a representation of the potential space using the reparameterization trick. The calculation formula is as follows:

[0081] For the first-layer graph convolution:

[0082] The l-th layer graph convolution is:

[0083] Calculate the mean and standard deviation of the potential space:

[0084]

[0085] μ is the mean, and logσ is the standard deviation.

[0086] Use the reparameterization trick:

[0087] Z = μ + ε·exp(logσ)

[0088] where ReLU is the ReLU activation function.

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

[0090] The present invention proposes the VGAE-CCI model, which can identify cell-cell communication in three-dimensional space. In the face of different levels of noise and missing spatial transcriptomics datasets, it can accurately identify cell-cell communication methods from missing or partially missing data, achieving high stability and accuracy. On multiple single-cell and spatial transcriptomics datasets, the VGAE-CCI model can achieve excellent cell-cell communication recognition. By comparing with other latest cell-cell communication recognition models, the VGAE-CCI model shows better performance in multiple evaluation metrics. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 is a schematic flowchart of a method for identifying cell communication based on multi-omics data according to the present invention;

[0092] Figure 2 is a structural model diagram of the VGAE-CCI model of the present invention;

[0093] Figure 3 is a schematic diagram of the comparison results between VGAE-CCI and other methods on the scRNA-seq dataset of the present invention, where AUROC is the area under the receiver operating characteristic curve;

[0094] Figure 4It is a schematic diagram of the performance of VGAE-CCI under different noises on scRNA-seq and ST-seq datasets of the present invention. Among them, seqFISH and MERFISH are two spatial transcriptomics datasets, std_dev is the standard deviation of the noise, epoch is the time point, and AP is the average precision;

[0095] Figure 5 It is a schematic diagram of the performance of VGAE-CCI of the present invention in reconstructing the cell interaction network. Among them, ACC is the accuracy;

[0096] Figure 6 It is a schematic diagram of the intercellular communication recognized by VGAE-CCI of the present invention across three layers of tissue sections 151670, 151671, and 151672 between tissues;

[0097] Figure 7 It is a schematic diagram of the visualization display of VGAE-CCI of the present invention on the dataset. Detailed implementation manners

[0098] Next, the present invention will be further described with reference to the accompanying drawings:

[0099] Detailed implementation manner 1: Combining Figure 1-7 To illustrate this implementation manner, as Figure 1 shown, the method for identifying cell communication based on multi-omics data described in this implementation manner is achieved through the following steps:

[0100] S1: Obtain single-cell and spatial transcriptomics datasets, and the single-cell and spatial transcriptomics datasets include datasets composed of scRNA-seq data and ST-seq data. Specifically, in this implementation manner, VGAE-CCI is compared with six methods for identifying cell communication, namely DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI, and DeepLinc, and ACC, AUROC, and AP metrics are used to evaluate the effects of these models in identifying cell communication on scRNA-seq and ST-seq datasets. We obtained 6 sets of real multi-omics datasets from the prior art. They are the mouse visual cortex dataset analyzed by seqFISH, the mouse hypothalamus slice dataset analyzed by MERFISH technology, the scRNA-seq dataset of human testis, the spatial transcriptomics dataset of human colorectal cancer liver metastasis, the human melanoma dataset, and the human dorsolateral prefrontal cortex dataset analyzed by the Visium platform.

[0101] The single-cell and spatial transcript data datasets obtained in this implementation manner are as follows:

[0102] The seqFISH dataset data was processed according to the procedures described in the original study. The final dataset contained 1,597 cells and 125 genes.

[0103] For the MERFISH dataset, we used the section area at Bregma +0.11 mm from animal No. 18 because this area contained the most single cells. After removing blurred cells, the final dataset contained 4,975 cells.

[0104] The scRNA-seq dataset of the human testis contains 3,074 cells.

[0105] The spatial transcriptomics dataset of human colorectal cancer (CRC) liver metastases contains four different cancer patients. Here, we selected patient 1, which contains 2,887 cells.

[0106] The human melanoma dataset contains 4,645 cells.

[0107] For the human dorsolateral prefrontal cortex (DLPFC) dataset analyzed by the Visium platform, we used three sections, 151670, 151671, and 151672, which contained a total of 11,465 cells.

[0108] S2: Preprocess the single-cell and spatial transcriptomics datasets described in S1. For scRNA-seq data, first perform data cleaning to remove cells and genes with all-zero values as a preliminary screening step. Subsequently, normalize the screened gene expression matrix to reduce data bias and improve analysis accuracy. For ST-seq data, we obtain the corresponding spatial coordinates and then construct the corresponding cell adjacency matrix.

[0109] Accurately predicting biologically meaningful cell-cell communication requires comprehensive signal L-R pairs. In recent years, many published L-R pair databases have emerged, such as CellChatDB, CellTalkDB, CellCallDB, CellPhoneDB, connectomeDB2020, and iTALK, etc. Their emergence is all for better predicting CCI, and these databases provide literature-supported L-R pairs for humans and mice.

[0110] To better predict cell-cell communication, we integrated publicly available ligand-receptor pair databases to construct the SCTDB database. In SCTDB, human and mouse L-R interactions are included. To ensure that the data in SCTDB is based on reliable and validated scientific research, we prioritized interaction pairs supported by the literature to improve the accuracy and reliability of predicting CCI. In addition, we also considered the existence of multi-subunit complexes as they are crucial functional units in cellular processes and involve complex L-R interactions. By considering these factors, we ensure that the predicted interactions reflect the complexity in actual biological processes. Finally, we integrated all the data and removed redundant L-R pairs to ensure the uniqueness and accuracy of the dataset. Therefore, SCTDB retains 5,786 human-verified L–R interactions and 4,806 mouse-verified L–R interactions.

[0111] The VGAE-CCI framework takes the scRNA-seq gene expression matrix and the corresponding adjacency matrix as the main inputs. For the gene expression matrix, we first perform data cleaning to remove cells and genes with all-zero values as a preliminary screening step. Subsequently, we normalize the screened gene expression matrix to reduce data bias and improve the analysis accuracy. The adjacency matrix can be obtained by two methods: one is to obtain cell coordinate data and refer to the construction of the adjacency matrix in the method section; the other is based on real receptor-ligand pairs. Through SCTDB, the processed gene expression matrix is matched with SCTDB, duplicate genes are extracted to construct a symmetric matrix, and finally, the adjacency matrix is constructed according to the LR pairs.

[0112] S3: Construct a deep neural network model. The deep neural network model includes an encoder, a variational inference module, a decoder, and an adversarial regularization module. The encoder is responsible for encoding the input graph data into a low-dimensional latent space. The variational inference module introduces a probabilistic approach to the embeddings generated by the encoder. The decoder is responsible for reconstructing the adjacency matrix of the graph or the probability of the existence of edges from the node embeddings in the latent space. The adversarial regularization module enhances the robustness of the model to input perturbations and noise.

[0113] Specifically, as Figure 2 shown, the VGAE-CCI model framework is used to identify cell communication in single-cell and spatial transcriptomic data. The encoder is a 3-layer graph convolutional network, and the decoder is a sigmoid function of the dot product of latent variables. First, an undirected adjacency network of cells is constructed using the SCTDB database or the spatial coordinates of cells, where nodes represent cells (or genes) and edges represent adjacent cell pairs (or L-R pairs), as Figure 2(a). The network is represented by the adjacency matrix A, and the single-cell gene expression data serves as the features of the nodes in the adjacency matrix. Subsequently, we input this adjacency network with node features into a neural network composed of three GCNs, as shown in Figure 2 (b). In the GCN, the generated latent space representation Z captures the feature information of individual cells and their neighboring cells, and Z is constrained by an adversarial regularization module from a prior Gaussian distribution. Based on this latent space representation, the decoder generates the adjacency matrix A' of the reconstructed intercellular interaction network through dot product operations, demonstrating the reconstructed intercellular interaction network. However, to achieve cross-tissue cell communication analysis, the ST-seq data and scRNA-seq data need to be aligned by the PASTE method. After alignment, we obtain three-dimensional spatial coordinates and three-dimensional spatial structure diagrams. Next, we input these three-dimensional coordinates, spatial structure diagrams, and gene expression matrices into the GCN to reconstruct the three-dimensional cell communication structure.

[0114] 1. S4: Train the deep neural network model constructed in S3 based on the preprocessed single-cell and spatial transcriptomics datasets in S2. The steps for training the deep neural network model include:

[0115] S401: Alignment of multi-layer tissue sections. To explore the analysis of cell communication between multi-layer tissue sections, the most important step is to align the multi-layer tissue sections to facilitate subsequent analysis and processing. Here, we use the mainstream PASTE algorithm to align the multi-layer tissue sections. This algorithm can better handle the noise and variability in the data compared to other algorithms, thereby increasing the robustness of the alignment, and it is more suitable for complex spatial transcriptomics data.

[0116] First, we extract feature vectors from each section, usually constructed using the gene expression data or coordinates of each cell. Next, calculate the distance between each cell. Here, the Euclidean distance is used to calculate the distance between cells. Assume xi and xj are the feature vectors of two cells, then their Euclidean distance dij is calculated as follows:

[0117]

[0118] where, x ik and x jk are the k-th components in the feature vector respectively.

[0119] Next, align the sections by optimizing the alignment objective function. First, initialize the matching matrix P:

[0120]

[0121] Among them, n and m are the numbers of points in the two slices respectively.

[0122] Then for each slice, calculate the distance matrix between its points:

[0123] D X (i,j) = ||X i - X j ||

[0124] D Y (k,l) = ||Y k - Y l ||

[0125] Among them, D X (i,j) and D Y (k,l) are the distance matrices of the two cell graphs respectively, and X and Y are the point sets of the two slices. The PASTE algorithm uses the Gromov-Wasserstein (GW) distance to align the two cell graphs. The GW distance is a measure of the difference between the distance matrices between cells, so as to minimize the GW distance between the two cell graphs. Its optimization objective function can be expressed as:

[0126] C(i,k,j,l) = (D X (i,j) - D Y (k,l)) 2

[0127]

[0128] Use the Sinkhorn-Knopp algorithm for iterative optimization to update the matching matrix P.

[0129] P ← diag(u)·P·diag(v)

[0130] Among them, u and v are vectors for normalization, and P is updated iteratively to meet the row and column sum constraint conditions.

[0131] Finally, use the generalized Procrustes analysis to optimize the rotation and scaling parameters to make the alignment result more accurate.

[0132]

[0133] Among them, R is the rotation matrix and t is the translation vector.

[0134] Through the above steps, the multi-layer tissue slices can be successfully aligned.

[0135] S402: Construction of the adjacency matrix, and the construction of the adjacency matrix is the construction of the communication adjacency matrix at the cell level or the construction of the communication adjacency matrix at the gene level.

[0136] When constructing the adjacency matrix, we provide two methods. When wanting to explore cell-level communication, coordinate data of cells need to be provided to construct the adjacency matrix; while when wanting to explore gene-level communication, the adjacency matrix can be constructed according to the provided real SCTDB database.

[0137] For the exploration of cell-level communication, we have the following method to construct the adjacency matrix (similar to DeepLinc). From a general biological perspective, it can be reasonably assumed that in solid tissues, most cells can directly contact three or more other cells, and generally, it is considered that the cell interacts with its neighboring cells. Therefore, we can construct the adjacency matrix based on this. First, we calculate the distance matrix according to the coordinate data of cells. Here, the Euclidean distance is used to calculate the distance between cells:

[0138]

[0139] where D (i,j) is the distance between cell i and cell j, and x ik and x jk are the coordinates of cell i and cell j in the k-th dimension.

[0140] Next, we determine the nearest neighbor cells of each cell. For each cell i, find its n nearest neighbor cells and store these neighbor cells in Ni:

[0141] N i ={j|D ij is the first n smallest elements in D i}

[0142] Here, D i represents the i-th row of the distance matrix D.

[0143] Finally, construct the adjacency matrix:

[0144]

[0145] Finally, we need to symmetrize the constructed adjacency matrix to make it a symmetric matrix.

[0146] When wanting to explore gene-level communication, the following method can be used to construct the adjacency matrix. If the given gene expression data is E = {e1, e2,..., eu}, where ei is the identifier of the i-th gene and u is the total number of genes. By combining the sorted LR pairs, an adjacency matrix A can be constructed. If genes i and j in the given gene expression matrix exist in the SCTDB and correspond to each other, then ei and ej are considered related, and the adjacency matrix Aij can be expressed as:

[0147]

[0148] Similarly, the constructed adjacency matrix should also be symmetrized at this time.

[0149] S403: Construction of the cell communication model. Combine GCN and VAE to construct a cell communication model VGAE-CCI based on VGAE.

[0150] Here, we construct a cell communication model based on VGAE. By combining GCN and VAE, an end-to-end model can be constructed to simultaneously learn the graph structure information and the latent representation from the cell communication network. The combination of the two can enhance the model's representation ability and generation ability for the cell communication network. And the cell communication network is usually sparse and high-dimensional, and the combination of GCN and VAE can effectively process such data and improve the robustness and generalization ability of the model.

[0151] VAE is a generative model that can learn the latent representation of data and generate new samples. This can be used in the cell communication network to generate potential receptor-ligand pairs and simulate new cell communication networks. In the encoder, when receiving the input data x and outputting the mean and standard deviation

[0152]

[0153] where and are learned through the neural network parameters , represents the normal distribution, represents the posterior distribution, which describes the distribution of the latent variable z given the input x.

[0154] To enable the gradient to be passed through the latent variable z, the reparameterization trick is used:

[0155] z = μ φ (x) + σ φ (x) ⊙ ∈

[0156] where is the noise sampled from the standard normal distribution.

[0157] The loss function includes the reconstruction loss and the KL divergence loss. The reconstruction loss is used to measure the difference between the generated data and the original data. Here, we use the cross-entropy loss:

[0158]

[0159] The KL divergence loss is used to measure the difference between the approximate posterior distribution and the prior distribution p(z):

[0160]

[0161] The total loss function is the weighted sum of the reconstruction loss and the KL divergence loss:

[0162]

[0163] where β is a weight parameter used to balance the reconstruction loss and the KL divergence loss. Here, the default setting is 1.

[0164] The GCN, as an encoder, adopts a method of aggregating features based on neighboring nodes, updates the target node by weighting the neighboring nodes, and the propagation method between layers in the GCN neural network layer is:

[0165]

[0166] In the formula: is the degree matrix, and the formula is H is the feature of each layer, and σ is the sigmoid function.

[0167] This propagation method can achieve parameter sharing for all nodes, that is, the update rules for each node are the same. In this way, all genes in the network are identified and classified according to the expression values, which is more in line with the needs.

[0168] In our model, features are extracted through the convolutional layer, and the representation in the latent space is generated through the reparameterization trick. For the first-layer graph convolution:

[0169]

[0170] The l-th layer graph convolution is:

[0171]

[0172] Next, calculate the mean and standard deviation of the latent space:

[0173]

[0174] μ is the mean, and logσ is the standard deviation.

[0175] For the reparameterization method:

[0176] Z = μ + ε·exp(logσ)

[0177] where

[0178] Finally, we use the inner product decoder to calculate the reconstructed adjacency matrix:

[0179]

[0180] where σ represents the sigmoid function.

[0181] S5: Identify cell-cell communication for the data to be tested based on the deep neural network model trained in S4. After obtaining single-cell and spatial transcriptome data and successfully constructing the model, we use VGAE-CCI for cell-cell communication identification.

[0182] This embodiment compares VGAE-CCI with six methods for identifying cell-cell communication, namely DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI, and DeepLinc. The accuracy (ACC), area under the receiver operating characteristic curve (AUROC), and average precision (AP) metrics are used to evaluate the performance of these models in identifying cell communication on scRNA-seq and ST-seq datasets. Among them, the seqFISH dataset was downloaded from the data portal, the MERFISH dataset was downloaded for the anterior hypothalamic area of the mouse, and we used the sliced area at Bregma +0.11 mm from animal No. 18. The scRNA-seq dataset of human testis (GSE106487). The spatial transcriptomics dataset of human colorectal cancer (CRC) liver metastasis contains four different cancer patients. Here, we selected patient 1 (GSE217414). The human melanoma dataset (GSE72056). The human dorsolateral prefrontal cortex (DLPFC) dataset can be downloaded from the portal.

[0183] Example 1

[0184] The spatial transcriptomics dataset of human colorectal cancer liver metastasis is a comprehensive dataset for studying high-throughput spatial gene expression during the liver metastasis of colorectal cancer. By analyzing cell-cell communication in this dataset, we can deeply understand the spatial distribution, migration path of cancer cells in the liver, and their interactions with other cell types in the microenvironment. This is of great significance for revealing the complex dynamics of the tumor microenvironment, the metastasis mechanism of cancer cells, and their potential therapeutic targets. Therefore, we compared the performance of VGAE-CCI with the other six methods on this dataset. Each method uses its default threshold, and the performance is evaluated by accuracy (ACC). VGAE-CCI achieved the highest ACC value, and the results are as follows:

[0185] Table 1 Comparison of ACC between VGAE-CCI and the other 6 methods on the CRC real dataset

[0186]

[0187] Example 2

[0188] It was compared with other methods on the human melanoma dataset. Similarly, we used ACC for performance evaluation, and VGAE-CCI also achieved the highest ACC value. The results are as follows:

[0189] Table 2 Comparison of ACC of VGAE-CCI with 6 other methods on the Melanoma real dataset

[0190]

[0191] The results of Example 1 and Example 2 show that VGAE-CCI can efficiently identify the corresponding L-R pairs in cancer datasets and exhibits high accuracy.

[0192] Example 3

[0193] As Figure 3 shown, the performance of VGAE-CCI was compared with that of 6 other CCI prediction methods on the human testicular cell dataset. Each method used its own default cut-off value. From Sertoli cells to SSCs, VGAE-CCI identified 69 cell-cell communications supported by the literature. While DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI, and DeepLinc identified 42, 42, 59, 40, 62, and 55 cell-cell communications supported by the literature, respectively. It was reported that among the cell-cell communications supported by the literature, 87.34% of the cell-cell communications identified by VGAE-CCI (69 / 79) were involved in spermatogenesis. The identification results of DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI, and DeepLinc are as Figure 3 shown in (a). Then we used AUROC to compare these methods, and as Figure 3 shown in the results of (b), VGAE-CCI achieved the highest area under the receiver operating characteristic (AUROC) curve. These results indicate that VGAE-CCI may infer cell-cell communications more accurately than these existing methods and has the ability to discover biologically significant CCI from scRNA-seq data.

[0194] Example 4

[0195] As Figure 4As shown, to demonstrate the performance of VGAE-CCI on noisy datasets, we applied VGAE-CCI to two spatial transcriptomics datasets, seqFISH and MERFISH. First, to test the tolerance of VGAE-CCI to artificial noise, specifically, we added Gaussian noise of different levels to the original gene expression data. On both datasets, it can be observed that under Gaussian noise conditions with a standard deviation between 0 and 7, the AUROC and AP of VGAE-CCI increased rapidly within the first 10 epochs (time points), and then the growth rate gradually became flat. In addition, regardless of the magnitude of the noise standard deviation, the values of AUROC and AP of VGAE-CCI ultimately showed high consistency. Therefore, VGAE-CCI has a high anti-interference ability, and AUROC and AP remain basically unchanged for different levels of noise (where std_dev represents the standard deviation of the noise). This indicates that VGAE-CCI is highly robust to noise in gene expression data, which is crucial when dealing with high-noise spatial transcriptome data.

[0196] Example 5

[0197] As Figure 5 shown, we evaluated the performance of VGAE-CCI in recovering a partially missing cell-cell interaction network. Specifically, we randomly removed different proportions of existing edges from the original cell-cell interaction network, and then used the remaining network and single-cell transcriptomics data to train the VGAE-CCI model. By calculating the values of AUROC, AP, and ACC, we evaluated the performance of VGAE-CCI in identifying arbitrarily removed different proportions of interactions. VGAE-CCI showed high precision in imputing missing interactions. Even when half of the interactions were missing in the input, VGAE-CCI was still able to recover these missing edges with an accuracy of 75-80%. This result indicates that VGAE-CCI can learn from a limited and incomplete intercellular interaction network and effectively recover the missing edges. This imputation and prediction ability of VGAE-CCI is crucial because it can infer the complete intercellular interaction landscape from missing or partially missing spatial transcriptome data.

[0198] The results of Examples 1 to 5 show that VGAE-CCI demonstrated high robustness and recovery ability for a partially missing cell-cell interaction network in gene expression data. Even when a large number of interactions were missing, it could still maintain high AUROC, AP, and ACC, proving its effectiveness in dealing with complex and incomplete spatial transcriptome data.

[0199] Example 6

[0200] VGAE-CCI identifies cell-cell communication in three-dimensional space. In the real physiological environment, cell-cell communication is usually three-dimensional. Single-layer tissue sections may miss the cell-cell communication information between layers during analysis. Therefore, it is particularly important to identify cell communication between multi-layer tissue sections. The complexity of this three-dimensional tissue structure is crucial for understanding the dynamic regulation of cell interactions and functions. Multi-layer sections can more accurately simulate the real interactions of cells in the whole tissue and reflect the behavior of cells in a complex physiological environment. The prerequisite for exploring cell-cell communication between multi-layer tissue sections is to align the multi-layer tissue sections.

[0201] Tissue sections usually have different spatial structures. Without alignment, the positions of cells on the sections will no longer be consistent with the actual biological tissue structure. This will lead to distortion of the relative positions and spatial relationships between cells. As a result, the spatial dependence of cell-cell communication may not be correctly reflected, which may affect the accuracy and interpretation of the cell communication pattern. As shown in Table 3, the effects before and after aligning the tissue sections are presented. Among them, Seed (random seed) is a fixed value used to control the generation of random numbers, which can ensure that the random behavior during the experiment and model training process is repeatable each time it runs. It can be seen that if the tissue sections are not aligned, the accuracy rate of the prediction results will decrease, which may lead to misunderstandings of cell functions, tissue structures, and biological processes, thus affecting research conclusions and clinical applications.

[0202] Table 3 Effects before and after aligning tissue sections

[0203]

[0204] As Figure 6 shown, we used three tissue sections 151670, 151671, and 151672 in the human dorsolateral prefrontal cortex dataset. As Figure 6 shown in (a) therein, the changes before and after alignment among the three-layer tissue sections are presented. To further understand the cell-cell communication between different sections, we applied the data to VGAE-CCI. Figure 6 In (b), (c), and (d) therein, the spatial position views of the cells after alignment and the views of cell-cell communication among the three layers of sections are shown. It can be seen that cell-cell communication does not only occur within a single section, but there are also cells that communicate across tissues with cells in other tissues.

[0205] Example 7

[0206] To further understand the behavior of interlayer communication, we visualized the communication network diagrams of the top 10 genes in different layers, as Figure 6(e) As shown. The ITGAV (Integrin alpha V) gene encodes an integrin protein and belongs to the integrin family. It is closely related to various cellular regulations, cell-matrix interactions, and biological processes such as the development of the skeletal system, immune response, and cell migration. Moreover, ITGAV is related to the occurrence and development of various diseases (especially cancer and autoimmune diseases). It plays an important role in tumor metastasis, inflammatory response, and platelet function. As Figure 5 in e, ITGAV is a gene shared by the three figures, but in different layers, the genes associated with it are different. Thus, it can be seen that VGAE-CCI can recognize these relationships, and can understand the relationships and influences between different layers, thereby exploring the intercellular communication mechanism between layers, which will help to better explore the differences between tumor tissues and normal tissues.

[0207] Example 8

[0208] VGAE-CCI can discover L-R pairs related to cancer. Colorectal cancer refers to a malignant tumor that occurs in the colon or rectum and has a high incidence and fatality rate globally. Liver metastasis of CRC (CRCLM) is one of the most common metastasis forms of CRC, which means that cancer cells spread from the primary tumor site in the colon or rectum through the blood or lymphatic system and invade the liver tissue. In this study, we used the spatial transcriptome data of colorectal cancer liver metastasis to predict and analyze intercellular communication. Using an advanced computational model, we inferred the ligand-receptor interactions between different cells and identified gene pairs with significant communication potential during liver metastasis. Figure 7Figures (a)-(d) show the network diagrams of the top 10 gene pairs with the highest communication probability during the analysis of colorectal cancer liver metastasis data by VGAE-CCI. Each node represents a gene, and each edge indicates a high communication probability between two genes. Through such network diagrams, we can visually observe the communication patterns of specific gene pairs in the liver metastasis microenvironment, thus providing important references for revealing the complex dynamics of the tumor microenvironment and discovering potential therapeutic targets. For example, EGFR (epidermal growth factor receptor) is a gene encoding the epidermal growth factor receptor, and its importance in cancer has been widely recognized. The role of EGFR in colorectal cancer liver metastasis is particularly crucial. Its activation can initiate multiple intracellular signaling pathways, such as the PI3K / AKT, RAS / RAF / MEK / ERK, and JAK / STAT pathways, and then affect biological processes such as cell proliferation, apoptosis, angiogenesis, and migration. PLD2 (phospholipase D2) is also involved in processes such as cell signaling, cell migration and invasion, and cell proliferation. PLD2 mediates intracellular signaling by generating the second messenger phosphatidic acid (PA), and interacts with multiple signaling molecules to regulate cell physiological activities. EGFR and PLD2 interact with each other in cell signaling and synergistically promote cell proliferation, migration, and invasion. In the EGFR signaling pathway, PLD2 is considered one of the downstream effector molecules, which enhances the cell's response to EGFR-mediated growth factors by regulating key nodes in the EGFR signaling pathway. Their synergistic effect plays a key role in tumor progression, and understanding this interaction is of great significance for developing new treatment strategies and overcoming drug resistance. Research shows that by targeting EGFR and PLD2 simultaneously, the drug resistance problem of single-target drug therapy can be overcome. Therefore, in-depth exploration of the interaction mechanism between EGFR and PLD2 in colorectal cancer liver metastasis is of great significance for tumor biology research and the formulation of clinical treatment strategies.

[0209] Example 9

[0210] Melanoma is a malignant tumor derived from skin melanocytes (i.e., pigment-producing cells). Although it accounts for only a small part of all skin cancers, it is the most lethal one, with high invasiveness and metastasis. In the human melanoma dataset, we visualized the L-R pairs of communication between different cells (such as Figure 7(as shown in (a)-(d)). In the onset, prognosis, and treatment of melanoma, it is closely related to gene mutations such as BRAF, CKIT, and NRAS. The overall incidence of NRAS mutations in melanoma is approximately 20-25%, and NRAS mutations participate in signal transduction through pathways such as RAF-MEK-ERK and P13K / Akt / mTOR, thereby affecting cell growth and differentiation as well as the occurrence and development of tumors. The BRAF gene encodes the serine / threonine protein kinase B-Raf, which is part of the MAPK / ERK signaling pathway and participates in the regulation of cell division, differentiation, and apoptosis. Among them, the BRAF V600E mutation is the most common melanoma-related mutation, which continuously activates the B-Raf kinase, leading to excessive cell proliferation and tumor formation. The CDKN2A gene encodes two proteins: p16INK4a and p14ARF, both of which are tumor suppressors that prevent cell cycle progression and promote apoptosis by inhibiting cyclin-dependent kinases (CDK4 / 6) and stabilizing p53, respectively. CDKN2A gene mutations or deletions will result in the loss of function of p16 and p14, thus losing control of the cell cycle and promoting cell proliferation and tumor formation. As Figure 7 (as shown in (e)), which shows the communication relationships between BRAF, CDKN2A, and NRAS and other genes. In melanoma, BRAF and NRAS mutations promote cell proliferation by activating the MAPK / ERK pathway, while CDKN2A mutations further enhance the proliferation ability of tumor cells by relieving the inhibition of the cell cycle. When these three genes are mutated simultaneously, tumor cells will gain a stronger growth and survival advantage, leading to the rapid development and deterioration of melanoma.

[0211] In summary, VGAE-CCI plays a key role in elucidating the complex intercellular interactions and communication networks in various cancer types (such as colorectal cancer liver metastasis and melanoma). By using advanced computational models, VGAE-CCI can predict and visualize ligand-receptor interactions and identify genes that play key roles in the tumor microenvironment. This enables us to gain a deeper understanding of cell behavior and interactions, which is crucial for discovering potential therapeutic targets and developing effective treatment strategies.

[0212] The above are only the preferred embodiments of the present invention, and do not impose any form of limitation on the present invention. Although the present invention has been disclosed above with the preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to equivalent embodiments by using the disclosed technical content within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention, any simple modification, equivalent replacement, and improvement made to the above embodiments within the spirit and principle of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for identifying cell communication based on multi-omics data, characterized in that: The method for identifying cell communication based on multi-omics data is implemented through the following steps: S1: Obtain single-cell and spatial transcriptomics datasets, where the single-cell and spatial transcriptomics datasets include a dataset composed of scRNA-seq data and ST-seq data; S2: Preprocess the single-cell and spatial transcriptomics datasets described in S1. For scRNA-seq data, first perform data cleaning to remove cells and genes with all-zero values as a preliminary screening step. Subsequently, standardize the screened gene expression matrix to reduce data bias and improve analysis accuracy. For ST-seq data, we obtain the corresponding spatial coordinates and then construct the corresponding cell adjacency matrix; S3: Construct a deep neural network model, which includes an encoder, a variational inference module, a decoder, and an adversarial regularization module. The encoder is responsible for encoding the input graph data into a low-dimensional latent space. The variational inference module introduces a probabilistic method for the embeddings generated by the encoder. The decoder is responsible for reconstructing the adjacency matrix of the graph or the probability of the existence of edges from the node embeddings in the latent space. The adversarial regularization module improves the robustness of the model to input perturbations and noise; S4: Train the deep neural network model constructed in S3 based on the preprocessed single-cell and spatial transcriptomics datasets in S2; The steps for training the deep neural network model described above include: S401: Alignment of multi-layer tissue sections; S402: Construction of the adjacency matrix, where the construction of the adjacency matrix is the construction of the cell-level communication adjacency matrix or the gene-level communication adjacency matrix; S403: Construction of the cell communication model, combine GCN and VAE to construct a cell communication model VGAE-CCI based on VGAE; The steps for constructing the cell communication model are as follows: a. Combine GCN and VAE to construct an end-to-end model, and simultaneously learn the graph structure information and the latent space representation from the cell communication network. From the output of the encoder, the model learns to generate the distribution of the latent variable z, and the expression is as follows: where x is the input data, represents a normal distribution, and are the mean and standard deviation learned through the neural network parameters respectively, and represents that the posterior distribution describes the distribution of the latent variable z given the input x; b. Introduce a probabilistic method and use the reparameterization trick to learn the mean and variance of the latent distribution, so that the gradient can be passed through the latent variable z, and the expression is as follows: z = μ φ (x) + σ φ (x) ⊙ ∈ wherein, is the noise sampled from the standard normal distribution; c. Use the inner product decoder to calculate the reconstructed adjacency matrix, and the expression is as follows: Among them, is the sum of the identity matrix and the relationship matrix of each node, σ is the sigmoid function, Z is the latent space representation, and Z is constrained by the adversarial regularization module from the prior Gaussian distribution, and Z T is the transpose of Z; S5: Identify intercellular communication for the data to be measured based on the deep neural network model trained in S4.

2. The method for identifying cell communication based on multi-omics data according to claim 1, wherein: The alignment method described in S401 is the PASTE algorithm, and the alignment method steps of the PASTE algorithm include: S401.1: Extract feature vectors from each slice, and construct them using the gene expression data or coordinates of each cell; S401.2: Calculate the distance between each cell; S401.3: Align the slices by optimizing the alignment objective function. Calculate the distance matrix between the points for each slice; S401.4: Optimize the rotation and scaling parameters to make the alignment result more accurate.

3. The method for identifying cell communication based on multi-omics data according to claim 2, wherein: The calculation steps for each slice in S401.3 include: a. Initialize the matching matrix P, defined as follows: where n and m are the numbers of points in the two slices respectively; b. For each slice, calculate the distance matrix between its points, with the expression as follows: D X (i,j) = ||X i -X j || D Y (k, l) = ||Y k - Y l || Among them, D X (i, j) and D Y (k, l) are the distance matrices of two cell graphs respectively, and X and Y are the point sets of two slices; c. Use the Sinkhorn-Knopp algorithm for iterative optimization to update the matching matrix P, with the expression as follows: P←diag(u)·P·diag(v) where u and v are vectors for normalization, and P is iteratively updated to satisfy the row and column sum constraint conditions.

4. The method for identifying cell communication based on multi-omics data according to claim 1, wherein: The construction steps of the cell-level communication adjacency matrix described in S402 are as follows: a. Calculate the distance matrix according to the coordinate data of the cells. The calculation method is the Euclidean distance, and its calculation formula is as follows: where D (i,j) is the distance between cell i and cell j, and x ik and x jk are the coordinates of cell i and cell j in the k-th dimension; b. Determine the nearest neighbor cells of each cell. For each cell i, find its n nearest neighbor cells and store these neighbor cells in N i where N i is expressed as follows: N i = {j | D ij is the first n smallest elements in D i} where D i represents the i-th row of the distance matrix D; c. Construct the adjacency matrix, with the expression as follows: d. Symmetrize the constructed adjacency matrix to make it a symmetric matrix.

5. The method for identifying cell communication based on multi-omics data according to claim 1, characterized in that: The construction steps of the gene-level communication adjacency matrix described in S402 are as follows: a. Given gene expression data E = {e1, e2,..., eu}, where ei is the identifier of the i-th gene and u is the total number of genes, b. Construct the adjacency matrix, with the expression as follows: where if genes i and j in the given gene expression matrix exist in the database and correspond to each other, then ei and ej are considered related; c. Symmetrize the constructed adjacency matrix to make it a symmetric matrix.

6. The method for identifying cell communication based on multi-omics data according to claim 1, wherein: The loss function of the cell communication model VGAE-CCI is composed of the reconstruction loss and the KL divergence loss combined together and is defined as follows: where β is a weight parameter used to balance the reconstruction loss and the KL divergence loss, and the default setting is 1.

7. The method for identifying cell communication based on multi-omics data according to claim 6, characterized in that: The propagation method between layers in the GCN neural network layer is: Among them, is 's degree matrix, and the formula is H is the feature of each layer, and W is the learnable weight matrix of each layer.

8. The method for identifying cell communication based on multi-omics data according to claim 7, wherein: The potential space representation method is to extract features through a convolutional layer and generate the representation of the potential space through the reparameterization trick. The calculation formula is as follows: For the first layer of graph convolution: The l-th layer graph convolution is as follows: Calculate the mean and standard deviation of the potential space: μ is the mean and logσ is the standard deviation, Use the reparameterization trick: Z = μ + ε·exp(logσ) Among them ReLU is the ReLU activation function.

Images