Method for identifying cell communication based on multi-omics data

Abstract

The invention provides a method for identifying cell communication based on multi-omics data. The method is realized by the following steps: step 1, acquiring a single cell and space transcriptome data set; 2, preprocessing the single cell and the spatial transcriptomics data set; and step 3, constructing a deep neural network model. And step 4, training the constructed deep neural network model by using the preprocessed single cell and the space transcriptomics data set. And 5, carrying out intercellular communication identification on to-be-detected data by using the trained deep neural network model. According to the method, an excellent inter-cell communication recognition effect can be achieved on multiple sets of single-cell and spatial transcriptomics data sets. Compared with other latest intercellular communication identification methods, the method shows more excellent performance on a plurality of evaluation indexes. In addition, under the condition of different degrees of noise and on a missing space transcriptome data set, high stability and accuracy are achieved.

Description

Technical Field

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

[0002] Intercellular communication plays a vital role in maintaining the normal function of organisms, regulating development and differentiation, and controlling immune responses. Understanding and studying cell communication is of key significance to 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 junction proteins. The life of multicellular organisms depends on the coordination of cell activities, which in turn depends on cell-cell interactions (CCI) between different cell types and tissues in the organism. Therefore, the study of cell function increasingly needs to consider the spatial environment of each cell.

[0003] In recent years, several computational strategies based on ligand-receptor (LR) 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 analysis of gene expression and spatial organization at single-cell resolution, revealing the complex interactions between different cell types and subtypes. Single-cell sequencing technology can accurately identify the transcriptome characteristics of individual cells, revealing cell heterogeneity and dynamic changes, while spatial transcriptome sequencing technology provides spatial location information of cells in their original tissue environment, preserving tissue structure and spatial relationships between cells. Combining these two technologies, researchers can construct a more comprehensive and accurate cell communication network, revealing how cells communicate through ligand-receptor interactions, extracellular matrix components, and other signaling pathways.

[0004] However, spatial transcriptome data often contain technical noise and errors, such as inaccurate and missing data due to variations in sample processing and sequencing. Therefore, a method is needed to accurately identify intercellular communication from missing or partially missing data to improve the accuracy and reliability of cell communication research. Current research methods usually only predict intercellular communication on a single tissue slice, ignoring the connection between multi-layer tissue slices. Summary of the invention

[0005] The present invention aims to solve the problem of identifying intercellular communication in three-dimensional space, and further proposes 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:

[0007] 1. A method for identifying cell communication based on multi-omics data, wherein the method for identifying cell communication based on multi-omics data is implemented by the following steps:

[0008] S1: Acquire single-cell and spatial transcriptome datasets, wherein the single-cell and spatial transcriptome datasets include datasets consisting 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, we first perform data cleaning to remove cells and genes with all zero values ​​as a preliminary screening step. Then, we 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, 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 embedding generated by the encoder. The decoder is responsible for reconstructing the adjacency matrix of the graph or the existence probability of the edge from the node embedding in the latent space. The adversarial regularization module improves 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: Based on the deep neural network model trained in S4, the intercellular communication of the test data is identified.

[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: constructing an adjacency matrix, wherein the construction of the adjacency matrix is ​​the construction of an adjacency matrix for communication between cells or the construction of an adjacency matrix for communication between genes;

[0016] S403: Construction of the cell communication model, combining GCN with VAE to construct a VGAE-CCI cell communication network model based on VGAE.

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

[0018] S401.1: Extract feature vectors from each slice and construct them using gene expression data or coordinates of each cell;

[0019] S401.2: Calculate the distance between each cell;

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

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

[0022] Furthermore, the method for calculating the distance between cells in S401.2 is the Euclidean distance. Assuming that xi and xj are the feature vectors of two cells, the Euclidean distance dij between them is calculated as follows:

[0023]

[0024] Among them, x ik and x jk are the kth components in the eigenvector respectively.

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

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

[0027]

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

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

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

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

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

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

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

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

[0036] Furthermore, the optimization 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] GW(C,P)=∑ i，j，k，l C(i,k,j,l)P(i,k)P(j,l)

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

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

[0041]

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

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

[0044] a. Calculate the distance matrix based on the coordinate data of the cell. The calculation method is Euclidean distance, and the calculation formula is as follows:

[0045]

[0046] Among them, D (i，j) is the distance between cell i and cell j, x ik and x jk are the coordinates of cell i and cell j in the kth dimension;

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

[0048] N i ={j|D ij Yes D i The first n smallest elements in

[0049] Where Di represents the i-th row of the distance matrix D;

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

[0051]

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

[0053] Furthermore, the steps for constructing the gene-level communication adjacency matrix 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 the adjacency matrix, the expression is as follows:

[0056]

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

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

[0059] Furthermore, the steps for constructing the cell communication model in S403 are as follows:

[0060] a. Combine GCN with VAE to build an end-to-end model that simultaneously learns 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, which is expressed as follows:

[0061]

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

[0063] b. The probabilistic method uses the reparameterization technique to learn the mean and variance of the potential distribution so that the gradient can be passed through the latent variable z. The expression is as follows:

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

[0065] in, is the noise sampled from a standard normal distribution;

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

[0067]

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

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

[0070]

[0071] Among them, β is a weight parameter used to balance the reconstruction loss and KL divergence loss, and is set to 1 by default;

[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 log p θ (x|z) is the logarithm of the likelihood function, which represents the probability of reconstructing the input x given the latent variable z.

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

[0076]

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

[0078]

[0079] in, yes The degree matrix of H is the node feature matrix of each layer, and W is the learnable weight matrix of each layer.

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

[0081] For the first layer of graph convolution:

[0082] The lth layer of graph convolution is:

[0083] Compute the mean and standard deviation of the latent space:

[0084]

[0085] Use the reparameterization trick:

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

[0087] in ReLU is the ReLU activation function.

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

[0089] The present invention proposes a VGAE-CCI model, which can identify intercellular communication in three-dimensional space. In the face of varying degrees of noise and missing spatial transcriptome data sets, it can accurately identify intercellular communication methods from missing or partially missing data, achieving higher stability and accuracy. The present invention uses the VGAE-CCI model to achieve excellent intercellular communication recognition on multiple single-cell and spatial transcriptome data sets. By comparing with other latest models for identifying intercellular communication, the VGAE-CCI model shows better performance on multiple evaluation indicators. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] Figure 1 It is a schematic flow chart of a method for identifying cell communication based on multi-omics data of the present invention;

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

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

[0093] Figure 4Schematic diagram of the performance of VGAE-CCI under different noises on scRNA-seq and ST-seq datasets of the present invention, wherein seqFISH and MERFISH are two spatial transcriptome datasets, std_dev is the standard deviation of the noise, epoch is the time point, and AP is the average precision;

[0094] Figure 5 Schematic diagram of the performance of VGAE-CCI of the present invention in reconstructing cell interaction networks, where ACC is the accuracy rate;

[0095] Figure 6 This is a schematic diagram of the intercellular communication among three tissue slices 151670, 151671, and 151672 by VGAE-CCI cross-tissue recognition of the present invention;

[0096] Figure 7 It is a schematic diagram of the visualization display of VGAE-CCI on the data set of the present invention. DETAILED DESCRIPTION

[0097] Below, the present invention is further described in conjunction with the accompanying drawings:

[0098] Specific implementation method 1: Combination Figure 1-7 To illustrate this embodiment, Figure 1 As shown, the method for identifying cell communication based on multi-omics data described in this embodiment is achieved by the following steps:

[0099] S1: Obtain single-cell and spatial transcriptome datasets, which include datasets composed of scRNA-seq data and ST-seq data. Specifically, this embodiment compares VGAE-CCI with six methods for identifying intercellular communication: DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI, and DeepLinc, and uses ACC, AUROC, and AP indicators to evaluate the effects of these models on identifying cell communication on scRNA-seq and ST-seq datasets. We have 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 transcriptome dataset of human colorectal cancer liver metastasis, the human melanoma dataset, and the human dorsolateral prefrontal cortex dataset analyzed by the Visium platform.

[0100] The single cell and spatial transcriptome data sets obtained in this embodiment are as follows:

[0101] The seqFISH dataset data were processed according to the procedures described in the original study, and the final dataset contained 1597 cells and 125 genes.

[0102] For the MERFISH dataset, we used the slice region at Bregma + 0.11 mm from animal 18 because this region contained the most single cells. After removing ambiguous cells, the final dataset contained 4975 cells.

[0103] The scRNA-seq dataset of human testis contains 3074 cells.

[0104] The human colorectal cancer (CRC) liver metastasis spatial transcriptome dataset contains four different cancer patients. Here, we selected patient 1, which contains 2887 cells.

[0105] The human melanoma dataset contains 4645 cells.

[0106] The human dorsolateral prefrontal cortex (DLPFC) dataset analyzed by the Visium platform, we used three slices 151670, 151671 and 151672, containing a total of 11465 cells.

[0107] S2: The single-cell and spatial transcriptomics datasets described in S1 were preprocessed. For the scRNA-seq data, data cleaning was first performed to remove cells and genes with all values ​​of 0 as a preliminary screening step. Subsequently, the screened gene expression matrix was standardized to reduce data bias and improve analysis accuracy. For the ST-seq data, we obtained the corresponding spatial coordinates and then constructed the corresponding cell adjacency matrix.

[0108] Accurate prediction of biologically meaningful intercellular communication requires comprehensive signal LR pairs. In recent years, many published LR pair databases have emerged, such as CellChatDB, CellTalkDB, CellCallDB, CellPhoneDB, connectomeDB2020, and iTALK. They all appeared to better predict CCI, and these databases provide literature-supported LR pairs for humans and mice.

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

[0110] The VGAE-CCI framework takes the scRNA-seq gene expression matrix and the corresponding adjacency matrix as the main input. 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 standardize the screened gene expression matrix to reduce data bias and improve 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, matching the processed gene expression matrix with SCTDB, extracting duplicate genes to construct a symmetric matrix, and finally constructing the adjacency matrix based on the LR pairs.

[0111] 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 embedding generated by the encoder. The decoder is responsible for reconstructing the adjacency matrix of the graph or the existence probability of the edge from the node embedding in the latent space. The adversarial regularization module improves the robustness of the model to input perturbations and noise.

[0112] Specifically, Figure 2 As shown in Figure 1, the VGAE-CCI model framework is used to identify cell communication for single-cell and spatial transcriptome data. The encoder is a 3-layer graph convolutional network, and the decoder is a sigmoid function of the dot product of the 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 LR pairs), as shown in Figure 1. Figure 2(a). The network is represented by the adjacency matrix A, and the single-cell gene expression data is used as the features of the nodes in the adjacency matrix. Subsequently, we input the adjacency network with node features into the neural network composed of three GCNs, as shown in Figure 2 (b) As shown. In GCN, the generated latent space representation Z captures the characteristic information of a single cell and its 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 an adjacency matrix A' of the reconstructed intercellular interaction network through a dot product operation, showing the reconstructed intercellular interaction network. However, in order to achieve cross-tissue cell communication analysis, the ST-seq data and scRNA-seq data need to be aligned using the PASTE method. After alignment, we obtain three-dimensional spatial coordinates and a three-dimensional spatial structure map. Next, we input these three-dimensional coordinates and spatial structure maps into GCN in combination with the gene expression matrix to reconstruct the three-dimensional cell communication structure.

[0113] 2. S4: training the deep neural network model constructed in S3 based on the single cell and spatial transcriptomics datasets preprocessed by S2. The steps of training the deep neural network model include:

[0114] S401: Alignment of multi-layer tissue sections. In order 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 method PASTE algorithm to align multi-layer tissue sections. Compared with other algorithms, this algorithm can better handle the noise and variability in the data, thereby increasing the robustness of the alignment, and this algorithm is more suitable for complex spatial transcriptomics data.

[0115] First, we extract feature vectors from each slice, which are usually constructed using gene expression data or coordinates of each cell. Next, we calculate the distance between each cell. Here, we use the Euclidean distance to calculate the distance between cells. Assuming xi and xj are the feature vectors of two cells, the Euclidean distance dij between them is calculated as follows:

[0116]

[0117] Among them, x ik and x jk are the kth components in the eigenvector respectively.

[0118] Next, the slices are aligned by optimizing the alignment objective function. First, the matching matrix P is initialized:

[0119]

[0120] Where n and m are the number of points in the two slices, respectively.

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

[0122] D X (i, j) = ∥X i -X j ∥

[0123] D Y (k, l) = ∥Y k -Y l ∥

[0124] Among them, D X (i, j) and D Y (k, l) are the distance matrices of the two cell maps, 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 maps. The GW distance is a measure of the difference between the distance matrices between cells, thereby minimizing the GW distance between the two cell maps. Its optimization objective function can be expressed as:

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

[0126] GW(C,P)=∑ i，j，k，l C(i,k,j,l)P(i,k)P(j,l)

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

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

[0129] Among them, u and v are vectors used for normalization, and P is iteratively updated to satisfy the row, column and constraint conditions.

[0130] Finally, generalized Procrustes analysis is used to optimize the rotation and scaling parameters to make the alignment results more accurate.

[0131]

[0132] Where R is the rotation matrix and t is the translation vector.

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

[0134] S402: constructing an adjacency matrix, wherein the construction of the adjacency matrix is ​​the construction of an adjacency matrix for communication between cells or the construction of an adjacency matrix for communication between genes.

[0135] We provide two methods for constructing the adjacency matrix. When you want to explore cell-level communication, you need to provide the cell coordinate data to construct the adjacency matrix; when you want to explore gene-level communication, you can construct the adjacency matrix based on the real SCTDB database we provide.

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

[0137]

[0138] Among them, D (i，j) is the distance between cell i and cell j, x ik and x jk are the coordinates of cells i and j in the kth dimension.

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

[0140] N i ={j|D ij Yes D i The first n smallest elements in

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

[0142] Finally, construct the adjacency matrix:

[0143]

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

[0145] When you want to explore the communication between gene levels, you can use the following method to construct an adjacency matrix. If you are 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. By combining the sorted LR pairs, you can construct an adjacency matrix A. If there are genes i and j in the given gene expression matrix that exist in SCTDB and correspond to each other, then ei and ej are considered to be related, and the adjacency matrix Aij can be expressed as:

[0146]

[0147] Similarly, the constructed adjacency matrix should also be symmetric at this time.

[0148] S403: Construction of the cell communication model, combining GCN with VAE to construct a VGAE-CCI cell communication network model based on VGAE.

[0149] Here, we built a cell communication network model based on VGAE. By combining GCN with VAE, an end-to-end model can be built to simultaneously learn graph structure information and potential representations from the cell communication network. The combination of the two can enhance the model's representation and generation capabilities for cell communication networks. In addition, cell communication networks are usually sparse and high-dimensional. The combination of GCN and VAE can effectively process this type of data and improve the robustness and generalization ability of the model.

[0150] VAE is a generative model that can learn the latent representation of data and generate new samples. This can be used to generate potential receptor-ligand pairs and simulate new cell communication networks in cell communication networks. The encoder receives input data x and outputs the mean of the latent variable z. and standard deviation

[0151]

[0152] in, and is the neural network parameter What I have learned, represents a normal distribution, It means that the posterior distribution describes the distribution of the latent variable z given the input x.

[0153] To allow the gradient to propagate through the latent variable z, we use the reparameterization trick:

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

[0155] in, is the noise sampled from a standard normal distribution.

[0156] The loss function includes reconstruction loss and 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:

[0157]

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

[0159]

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

[0161]

[0162] Among them, β is a weight parameter used to balance the reconstruction loss and KL divergence loss. Here, it is set to 1 by default.

[0163] As an encoder, GCN uses a feature aggregation based on neighbor nodes, updates the target node by weighting the neighbor nodes, and the propagation method between layers in the GCN neural network layer is:

[0164]

[0165] In the formula: It is the sum of the identity matrix and the relationship matrix of each node. yes The degree matrix of H is the feature of each layer, and σ is the sigmoid function.

[0166] This propagation method can achieve parameter sharing among all nodes, that is, the update rules of each node are the same, so that all genes in the network can be identified and classified according to their expression values, which is more in line with the needs.

[0167] In our model, features are extracted through convolutional layers, and representations of the latent space are generated through reparameterization techniques. For the first layer of graph convolution:

[0168]

[0169] The l-th layer of graph convolution is:

[0170]

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

[0172]

[0173] For the reparameterized method:

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

[0175] in

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

[0177]

[0178] Where σ represents the sigmoid function.

[0179] S5: Based on the deep neural network model trained in S4, we identify intercellular communication on the test data. After obtaining single-cell and spatial transcriptome data and successfully building the model, we use VGAE-CCI to identify intercellular communication.

[0180] This embodiment compares VGAE-CCI with DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI, and DeepLinc, six methods for identifying intercellular communication, and uses the accuracy (ACC), area under the receiver operating characteristic curve (AUROC), and average precision (AP) indicators to evaluate the recognition of cell communication effects of these models on scRNA-seq and ST-seq datasets. Among them, the seqFISH dataset was downloaded from the data portal at: https: / / bitbucket.org / qzhu / smfish-hmrf / src / master / . The MERFISH dataset was downloaded from https: / / datadryad.org / stash / dataset / 10.5061 / dryad.8t8s24818 / . The section area at Bregma+0.11 mm from animal number 18 was used. scRNA-seq dataset of human testis (GSE106487). The human colorectal cancer (CRC) liver metastasis spatial transcriptome dataset contains four different cancer patients, here, we selected patient 1 (GSE217414). Human melanoma dataset (GSE72056). The human dorsolateral prefrontal cortex (DLPFC) dataset analyzed by the Visium platform can be downloaded from the portal at: http: / / spatial.libd.org / spatialLIBD / .

[0181] Example 1

[0182] The Human Colorectal Cancer Liver Metastasis Spatial Transcriptome Dataset is a comprehensive dataset for studying high-throughput spatial gene expression during colorectal cancer metastasis in the liver. By analyzing the intercellular communication in this dataset, we can gain a deep understanding of the spatial distribution of cancer cells in the liver, their migration pathways, 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 metastatic mechanisms of cancer cells, and their potential therapeutic targets. Therefore, we compared the performance of VGAE-CCI with six other methods on this dataset. Each method uses its default threshold and evaluates its performance by accuracy (ACC). VGAE-CCI achieved the highest ACC value, with the following results:

[0183] Table 1 Comparison of VGAE-CCI with ACC of other 6 methods on CRC real dataset

[0184]

[0185] Example 2

[0186] We compared it 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:

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

[0188]

[0189] The results of Example 1 and Example 2 show that VGAE-CCI can efficiently identify corresponding LR pairs in cancer datasets and exhibit a high accuracy rate.

[0190] Example 3

[0191] like Figure 3As shown in the figure, the performance of VGAE-CCI on the human testicular cell dataset is compared with that of 6 other CCI prediction methods. Each method uses its own default cutoff value. From Sertoli cells to SSCs, VGAE-CCI identified 69 literature-supported intercellular communications. DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI and DeepLinc identified 42, 42, 59, 40, 62 and 55 literature-supported CCIs, respectively. It is reported that among the literature-supported CCIs, 87.34% of the intercellular communications identified by VGAE-CCI (69 / 79) are involved in spermatogenesis. The identification results of DeepCCI, CellPhoneDB, SingleCellSignalR, CellChat, NATMI and DeepLinc are shown in the figure. Figure 3 Then we use AUROC to compare these methods, as shown in (a). Figure 3 In the results shown in (b), VGAE-CCI achieved the highest area under curve of receiver operating characteristic (AUROC). These results suggest that VGAE-CCI may be more accurate in inferring intercellular communication than these existing methods and has the ability to discover biologically meaningful CCIs from scRNA-seq data.

[0192] Example 4

[0193] like Figure 4 As shown in the figure, to demonstrate the performance of VGAE-CCI on noisy datasets, we applied VGAE-CCI to two spatial transcriptome datasets, seqFISH and MERFISH. First, in order to test the tolerance of VGAE-CCI to artificial noise, specifically, we added different degrees of Gaussian noise to the original gene expression data. On both datasets, it can be observed that under the condition of Gaussian noise with a standard deviation between 0 and 7, the AUROC and AP of VGAE-CCI rose rapidly within the first 10 epochs (time points), and then the growth rate gradually flattened. In addition, regardless of the size of the noise standard deviation, the AUROC and AP values ​​of VGAE-CCI finally showed high consistency. Therefore, VGAE-CCI has a high anti-interference ability, and AUROC and AP remain basically unchanged for different degrees of noise (where std_dev represents the standard deviation of the noise). This shows that VGAE-CCI has strong robustness to noise in gene expression data, which is crucial when processing high-noise spatial transcriptome data.

[0194] Example 5

[0195] like Figure 5 As shown, the performance of VGAE-CCI in recovering partially missing cell-cell interaction networks is evaluated. Specifically, we randomly remove different proportions of existing edges from the original cell-cell interaction network and then use the remaining network and single-cell transcriptome 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 interactions with different proportions removed arbitrarily. VGAE-CCI shows high accuracy in interpolating missing interactions. Even when half of the interactions are missing in the input, VGAE-CCI is still able to recover these missing edges with an accuracy of 75-80%. This result shows that VGAE-CCI is able to learn from limited and incomplete cell-cell interaction networks and effectively recover missing edges. This interpolation and prediction ability of VGAE-CCI is crucial because it can infer a complete cell-cell interaction landscape from missing or partially missing spatial transcriptome data.

[0196] The results of Examples 1 to 5 show that VGAE-CCI demonstrates high robustness and recovery capabilities for partially missing cell-cell interaction networks in gene expression data, and can maintain high AUROC, AP, and ACC even when a large number of interactions are lost, demonstrating its effectiveness in processing complex and incomplete spatial transcriptome data.

[0197] Example 6

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

[0199] Tissue slices usually have different spatial structures. If they are not aligned, the cell positions on the slices 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 communication may not be correctly reflected, which may affect the accuracy and interpretation of the communication pattern between cells. As shown in Table 3, the effects before and after aligning the tissue slices are shown. Among them, Seed is a fixed value used to control the generation of random numbers, which can ensure that the random behavior in the experiment and model training process is repeatable each time it is run. It can be seen that if the tissue slices are not aligned, the accuracy of the prediction results will decrease, which may lead to misunderstandings of cell functions, tissue structures and biological processes, thereby affecting research conclusions and clinical applications.

[0200] Table 3 Effects of aligning tissue sections before and after

[0201]

[0202] like Figure 6 As shown in Figure 2, we used three tissue sections 151670, 151671, and 151672 from the human dorsolateral prefrontal cortex dataset. Figure 6 As shown in (a), the changes before and after alignment of the three-layer tissue slices are shown. In order to further understand the communication between cells in different slices, we applied the data to VGAE-CCI. Figure 6 (b), (c) and (d) show the spatial position of aligned cells and the intercellular communication between three slices. It can be seen that intercellular communication is not just communication within a single slice, but also that cells can communicate across tissues with cells in other tissues.

[0203] Example 7

[0204] To further understand the behavior of inter-layer communication, we visualized the communication network diagram of the top 10 genes in different layers, such as Figure 6 (e) As shown. The ITGAV (integrin αv) gene encodes an integrin protein that belongs to the integrin family. It is closely related to various cell regulations, cell-matrix interactions, and biological processes such as the development of the skeletal system, immune response, and cell migration. Moreover, ITGAV is associated with the occurrence and development of many diseases (especially cancer and autoimmune diseases). It plays an important role in tumor metastasis, inflammatory response, and platelet function. Figure 5In the figure, ITGAV is a common gene in the three images, but the genes associated with it are different in different layers. It can be seen that VGAE-CCI can identify these relationships and understand the relationship and influence 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.

[0205] Example 8

[0206] VGAE-CCI can discover LR pairs associated with cancer. Colorectal cancer refers to malignant tumors that occur in the colon or rectum, with high incidence and mortality rates worldwide. Liver metastasis of CRC (CRCLM) is one of the most common forms of metastasis of CRC, which refers to the spread of cancer cells from the primary tumor site in the colon or rectum through the blood or lymphatic system and invasion of liver tissue. In this study, we used spatial transcriptome data of colorectal cancer liver metastasis to predict and analyze cell-to-cell communication. Using advanced computational models, we inferred ligand-receptor interactions between different cells and identified gene pairs with significant communication potential during liver metastasis. Figure 7(a)-(d) show the network graphs of the top 10 gene pairs with the highest communication probability in the process of VGAE-CCI analyzing colorectal cancer liver metastasis data. Each node represents a gene, and each edge represents a high communication probability between two genes. Through this network graph, we can intuitively observe the communication pattern of specific gene pairs in the liver metastasis microenvironment, thus providing an important reference 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 epidermal growth factor receptor, and its importance in cancer has been widely recognized. The role of EGFR in colorectal cancer liver metastasis is particularly critical. Its activation can initiate multiple intracellular signaling pathways, such as PI3K / AKT, RAS / RAF / MEK / ERK and JAK / STAT pathways, thereby affecting biological processes such as cell proliferation, apoptosis, angiogenesis and migration. PLD2 (phospholipase D2) is also involved in cell signaling, cell migration and invasion, and cell proliferation. PLD2 mediates intracellular signal transduction by generating the second messenger phosphatidic acid (PA), and interacts with a variety of signaling molecules to regulate cell physiological activities. EGFR and PLD2 interact in cell signal transduction and synergistically promote cell proliferation, migration and invasion. In the EGFR signaling pathway, PLD2 is considered to be one of the downstream effector molecules, which enhances the response of cells 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. Studies have shown that by simultaneously targeting EGFR and PLD2, 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.

[0207] Example 9

[0208] Melanoma is a malignant tumor that originates from melanocytes (pigment-producing cells) in the skin. Although it accounts for only a small proportion of all skin cancers, it is the most deadly one, being highly invasive and metastatic. In the human melanoma dataset, we visualized the LR pairs that communicate between different cells (e.g. Figure 7(a)-(d)). In the pathogenesis, prognosis and treatment of melanoma, BRAF, CKIT, NRAS and other gene mutations are closely related. The overall incidence of NRAS mutations in melanoma is about 20-25%. NRAS mutations participate in signal transduction through pathways such as RAF-MEK-ERK, 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, BRAF V600E mutation is the most common melanoma-related mutation, which causes the B-Raf kinase to be continuously activated, leading to excessive cell proliferation and tumor formation. The CDKN2A gene encodes two proteins: p16INK4a and p14ARF. Both proteins are tumor suppressors that prevent cell cycle progression and promote cell apoptosis by inhibiting cyclin-dependent kinases (CDK4 / 6) and stabilizing p53, respectively. CDKN2A gene mutation or deletion can lead to loss of function of p16 and p14, thereby losing control of the cell cycle and promoting cell proliferation and tumor formation. Figure 7 As shown in (e), the communication relationship between BRAF, CDKN2A and NRAS and other genes is demonstrated. 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 inhibitory effect on the cell cycle. When these three genes mutate at the same time, tumor cells will gain a stronger growth and survival advantage, leading to the rapid development and deterioration of melanoma.

[0209] In conclusion, VGAE-CCI has played a key role in elucidating the complex interactions and communication networks between cells in various cancer types, such as colorectal cancer liver metastasis and melanoma. By utilizing advanced computational models, VGAE-CCI is able to predict and visualize ligand-receptor interactions and identify genes that play key roles in the tumor microenvironment. This enables a deeper understanding of cell behavior and interactions, which is essential for discovering potential therapeutic targets and developing effective treatment strategies.

[0210] The above is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not used to limit the present invention. Any technician familiar with this profession can make some changes or modify the technical contents disclosed above into equivalent embodiments without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent replacement and improvement made to the above embodiments without departing from the content of the technical solution of the present invention, based on the technical essence of the present invention, within the spirit and principles 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 achieved by the following steps: S1: Acquire single-cell and spatial transcriptome datasets, wherein the single-cell and spatial transcriptome datasets include datasets consisting of scRNA-seq data and ST-seq data; S2: Preprocess the single-cell and spatial transcriptomics datasets described in S1. For scRNA-seq data, we first perform data cleaning to remove cells and genes with all zero values ​​as a preliminary screening step. Then, we 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 embedding generated by the encoder. The decoder is responsible for reconstructing the adjacency matrix of the graph or the existence probability of the edge from the node embedding 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 single-cell and spatial transcriptomics datasets preprocessed in S2; S5: Based on the deep neural network model trained in S4, the intercellular communication of the test data is identified.

2. A method for identifying cell communication based on multi-omics data according to claim 1, characterized in that: The steps of deep neural network model training described in S4 include: S401: Alignment of multi-layer tissue sections; S402: constructing an adjacency matrix, wherein the construction of the adjacency matrix is ​​the construction of an adjacency matrix for communication between cells or the construction of an adjacency matrix for communication between genes; S403: Construction of the cell communication model, combining GCN with VAE to construct a VGAE-CCI cell communication network model based on VGAE.

3. A method for identifying cell communication based on multi-omics data according to claim 2, characterized in that: 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 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, and for each slice, calculate the distance matrix between its points; S401.4: Optimize rotation and scaling parameters to make alignment results more accurate.

4. A method for identifying cell communication based on multi-omics data according to claim 3, characterized in that: Each slice calculation step in S401.3 includes: a. Initialize the matching matrix P, which is defined as follows: Where n and m are the number of points in the two slices respectively; b. For each slice, calculate the distance matrix between its points, expressed 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 the two cell maps, respectively, and X and Y are the point sets of the two slices; c. Use the Sinkhorn-Knopp algorithm for iterative optimization and update the matching matrix P. The expression is as follows: P←diag(u)·P·diag(v) Among them, u and v are vectors used for normalization, and P is iteratively updated to satisfy the row, column and constraint conditions.

5. The method for identifying cell communication based on multi-omics data according to claim 2, characterized in that: The steps for constructing the cell-level communication adjacency matrix in S402 are as follows: a. Calculate the distance matrix based on the coordinate data of the cell. The calculation method is Euclidean distance, and the calculation formula is as follows: Among them, D (i,j) is the distance between cell i and cell j, x ik and x jk are the coordinates of cell i and cell j in the kth dimension; b. Determine the nearest neighbor cells of each cell. For each cell i, find the n nearest neighbor cells and store these neighbor cells in N i In, N i The expression is as follows: N i ={j|D ij Yes D i The first n smallest elements in Where Di represents the i-th row of the distance matrix D; c. Construct the adjacency matrix, the expression is as follows: d. Symmetrize the constructed adjacency matrix to make it a symmetric matrix.

6. A method for identifying cell communication based on multi-omics data according to claim 2, characterized in that: The steps for constructing the gene-level communication adjacency matrix 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, the expression is as follows: Among them, if genes i and j in a given gene expression matrix exist in the database and correspond to each other, then ei and ej are considered to be related; c. Symmetrize the constructed adjacency matrix to make it a symmetric matrix.

7. A method for identifying cell communication based on multi-omics data according to claim 2, characterized in that: The steps for constructing the cell communication model described in S403 are as follows: a. Combine GCN with VAE to build an end-to-end model that simultaneously learns 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, which is expressed as follows: Among them, x is the input data, represents a normal distribution, and is the neural network parameter The learned mean and standard deviation, It means that the posterior distribution describes the distribution of the latent variable z given the input x; b. The probabilistic method uses the reparameterization technique to learn the mean and variance of the potential distribution so that the gradient can be passed through the latent variable z. The expression is as follows: z=μ φ (x)+σ φ (x)☉∈ in, is the noise sampled from a standard normal distribution; c. Use the inner product decoder to calculate the reconstructed adjacency matrix, the expression is as follows: in, is the sum of the identity matrix and the relationship matrix of each node, σ is the sigmoid function, Z is the latent space representation, Z is constrained by the adversarial regularization module from the prior Gaussian distribution, and Z T is the transpose of Z.

8. A method for identifying cell communication based on multi-omics data according to claim 7, characterized in that: The loss function of the cell communication network model VGAE-CCI By reconstruction loss and KL divergence loss Combined, defined as follows: Among them, β is a weight parameter used to balance the reconstruction loss and KL divergence loss, and is set to 1 by default.

9. A method for identifying cell communication based on multi-omics data according to claim 8, characterized in that: The propagation method between layers in the GCN neural network layer is: in, yes The degree matrix of H is the feature of each layer, and W is the learnable weight matrix of each layer.

10. A method for identifying cell communication based on multi-omics data according to claim 9, characterized in that: The latent space representation method extracts features through a convolutional layer and generates a representation of the latent space through a reparameterization technique. The calculation formula is as follows: For the first layer of graph convolution: The l-th layer of graph convolution is: Compute the mean and standard deviation of the latent space: Use the reparameterization trick: Z=μ+ε·exp(logσ) in ReLU is the ReLU activation function.

Images