Gene regulatory network inference method based on spatiotemporal transcriptomic data

Abstract

The present disclosure relates to the field of biological information, and in particular to a gene regulatory network inference method based on spatiotemporal transcriptomic data, and the use thereof. In the present disclosure, on the basis of a graph neural network, a graph variational auto-encoder is used, an adjacency matrix A is used as node connection information, cells and genes are used as node feature information, and position information of the cells is fully used, so as to reconstruct a feature matrix of nodes by means of three hidden layers, thereby performing more accurate cell clustering and dimensionality reduction analysis while inferring a gene regulatory network.

Description

A gene regulatory network inference method based on spatiotemporal transcriptome data Technical Field

[0001] The present disclosure relates to the field of bioinformatics, and in particular, to a gene regulatory network inference method based on spatiotemporal transcriptome data. Background Art

[0002] Transcription factors bind to genes through promoters, activating or inhibiting the subsequent transcription and translation of genes, thereby regulating gene expression. This regulatory relationship constitutes a gene regulatory network, which is of great significance for understanding the life processes of cells.

[0003] Summary of the Invention

[0004] With the development of spatiotemporal transcriptomics, compared to the gene expression matrix obtained from single-cell data, spatiotemporal transcriptomics can also provide cell location information. There are several gene regulatory network inference methods and algorithms:

[0005] (1) Based on pure mathematical methods. This type of method mainly uses pure mathematical modeling based on the gene regulatory network, and then solves the corresponding mathematical problems for inference. Representative methods include the Boolean model method (binarizing the gene expression matrix and inferring the gene regulatory network by optimizing the Boolean state of the model, such as BTR, etc.) and the differential equation method (obtaining pseudo time series based on embedding methods or other methods, and then constructing differential equations of the relationship between genes based on the time series, and obtaining the gene regulatory network by solving the equation, such as SCODE, etc.).

[0006] (2) Based on machine learning algorithms. This type of method mainly uses machine learning methods for inference, such as SCENIC. This method is based on a random forest regression model. It scores the relationship between transcription factors and target genes through the gene expression matrix, filters out modules with low scores, and thus constructs a gene regulatory network.

[0007] (3) Based on deep learning methods. This category mainly uses neural network construction for inference (such as DeepSEM). This method embeds the structural equation model into the variational autoencoder for inference.

[0008] Existing gene regulatory network inference methods are developed based on single-cell data, and there is currently a lack of methods for inferring gene regulatory networks based on spatiotemporal transcriptome data. The main problems with existing methods when applied to spatiotemporal transcriptome data are that they fail to fully utilize spatiotemporal information, leading to information loss during gene regulatory network inference, which reduces the accuracy of inference.

[0009] Therefore, there is an urgent need for a method that can infer gene regulatory networks based on spatiotemporal transcriptome data, especially a method that combines cell location information to infer gene regulatory networks.

[0010] The present disclosure aims to solve one of the technical problems in the related art at least to a certain extent.

[0011] To this end, the present disclosure provides, on one hand, a method for inferring gene regulatory networks based on spatiotemporal transcriptome data. According to an embodiment of the present disclosure, the method includes:

[0012] (A) Gene expression matrix and cell location matrix of samples were obtained using spatiotemporal transcriptomics platform;

[0013] (B) Using cells as feature dimensions and genes and cells as node dimensions, the gene expression matrix and the cell position matrix are merged to obtain a node feature matrix X with genes and cells as nodes and cells as features;

[0014] (C) Initialize the adjacency matrix A that represents the connection information of the graph nodes;

[0015] (D) inputting the initialized adjacency matrix A and node feature matrix X into the encoder part of the autoencoder, performing dimensionality reduction processing by the encoder, and extracting the embedding vector;

[0016] (E) inputting the embedding vector into the decoder part of the autoencoder to obtain a reconstructed node feature matrix Y;

[0017] (F) Based on the node feature matrix X and the reconstructed node feature matrix Y, obtain the trained adjacency matrix A1.

[0018] According to an embodiment of the present disclosure, the spatiotemporal transcriptome platform includes at least one selected from: Biotron S1000 spatial transcriptome sequencing platform, STOmics spatial transcriptome sequencing platform, and 10xgenomics spatial transcriptome sequencing platform.

[0019] According to an embodiment of the present disclosure, the sample includes a tissue sample selected from animals and plants.

[0020] According to an embodiment of the present disclosure, the adjacency matrix A includes three parts: gene-gene adjacency relationship, cell-cell adjacency relationship, and gene-cell adjacency relationship.

[0021] According to the embodiments of the present disclosure, a neural network uses the adjacency matrix A as a trainable parameter. Through learning and reconstruction by a graph variational autoencoder, the optimal control relationship is found. Because the graph variational autoencoder trains the adjacency matrix A, it needs to be initialized. According to the general process of neural network training, the adjacency matrix A is randomly initialized.

[0022] According to an embodiment of the present disclosure, in step (B), the direction of matrix merging is vertical merging.

[0023] According to an embodiment of the present disclosure, the autoencoder includes a graph variational autoencoder.

[0024] According to an embodiment of the present disclosure, the encoder part of the graph variational autoencoder uses a graph convolution layer for feature aggregation, and the decoder part of the graph variational autoencoder uses a graph convolution layer for reconstruction.

[0025] According to an embodiment of the present disclosure, in step (D), when the adjacency matrix A and the node feature matrix X are input into the encoder part, the hyperparameter n is introduced. channel The node feature matrix X of the two-dimensional tensor is expanded to a three-dimensional tensor, and after dimensionality reduction processing by the encoder, the embedding vectors of genes and cells are extracted respectively.

[0026] According to an embodiment of the present disclosure, in step (E), the reconstructed node feature matrix Y is obtained by the following method:

[0027] The embedding vector is input into the decoder of the graph variational autoencoder. By optimizing the loss function, the re-parameterization and Gumbel Softmax techniques are introduced to reconstruct the feature matrix Y of the node.

[0028] According to an embodiment of the present disclosure, the optimization of the loss function includes introducing an L1 regularization term for the adjacency matrix. Furthermore, the optimization of the loss function includes optimization in the process of reconstructing the feature matrix Y and optimization in the process of making the embedded vector obey a Gaussian distribution.

[0029] According to an embodiment of the present disclosure, after step (F), the method further includes step (G) of obtaining gene regulatory relationships based on the trained adjacency matrix A1.

[0030] According to an embodiment of the present disclosure, the adjacency matrix A, which is a trainable parameter, is optimized or updated based on the matrix Y and the matrix X (as training labels). The inter-gene relationships in the trained adjacency matrix A1 are the corresponding gene regulatory relationships.

[0031] In another aspect, the present disclosure provides a system for inferring gene regulatory networks based on spatiotemporal transcriptome data, comprising:

[0032] A data acquisition module, which is used to acquire the gene expression matrix and cell location matrix of the sample using the spatiotemporal transcriptome platform;

[0033] A matrix merging module, connected to the data acquisition module, uses cells as feature dimensions and genes and cells as node dimensions, and vertically merges the gene expression matrix and the cell position matrix to obtain a node feature matrix X with genes and cells as nodes and cells as features;

[0034] An initialization module, connected to the matrix merging module, for initializing an adjacency matrix A representing connection information of graph nodes;

[0035] A dimensionality reduction processing module, connected to the initialization module, is used to input the constructed adjacency matrix A and node feature matrix X into the encoder part of the graph variational autoencoder, and extract the embedding vector after the encoder dimensionality reduction processing;

[0036] A reconstruction module, connected to the dimensionality reduction processing module, configured to input the embedding vector into the decoder portion of a graph variational autoencoder to reconstruct a node feature matrix Y;

[0037] A feature training module is connected to the reconstruction module, and is used to obtain a trained adjacency matrix A1 based on the node feature matrix X and the reconstructed node feature matrix Y.

[0038] According to an embodiment of the present disclosure, in the matrix merging module, the direction of matrix merging is vertical merging.

[0039] According to an embodiment of the present disclosure, the spatiotemporal transcriptome platform includes at least one selected from the group consisting of the Biotron S1000 spatial transcriptome sequencing platform, the STOmics spatial transcriptome sequencing platform, and the 10xgenomics spatial transcriptome sequencing platform.

[0040] According to an embodiment of the present disclosure, the sample includes a tissue sample selected from animals and plants.

[0041] According to an embodiment of the present disclosure, the adjacency matrix A includes three parts: gene-gene adjacency relationship, cell-cell adjacency relationship, and gene-cell adjacency relationship.

[0042] According to an embodiment of the present disclosure, the autoencoder includes a graph variational autoencoder.

[0043] According to an embodiment of the present disclosure, the encoder part of the graph variational autoencoder uses a graph convolution layer for feature aggregation, and the decoder part of the graph variational autoencoder uses a graph convolution layer for reconstruction.

[0044] According to an embodiment of the present disclosure, in the dimensionality reduction processing module, when the adjacency matrix A and the node feature matrix X are input into the encoder part, the hyperparameter n is introduced. channel The node feature matrix X of the two-dimensional tensor is expanded to a three-dimensional tensor, and after dimensionality reduction processing by the encoder, the embedding vectors of genes and cells are extracted respectively.

[0045] According to an embodiment of the present disclosure, in the reconstruction module, the reconstructed node feature matrix Y is obtained by the following method:

[0046] The embedding vector is input into the decoder of the graph variational autoencoder. By optimizing the loss function, the re-parameterization and Gumbel Softmax techniques are introduced to reconstruct the feature matrix Y of the node.

[0047] According to an embodiment of the present disclosure, the optimization of the loss function includes introducing an L1 regularization term for the adjacency matrix. Furthermore, the optimization of the loss function includes optimization in the process of reconstructing the feature matrix Y and optimization in the process of making the embedded vector obey a Gaussian distribution.

[0048] According to an embodiment of the present disclosure, a data extraction module is further included after the feature training module. The data extraction module is connected to the feature training module, and the data extraction module obtains the gene regulatory relationship based on the trained adjacency matrix A1.

[0049] In another aspect, the present disclosure provides an electronic device, including a memory and a processor;

[0050] The processor runs a program corresponding to the executable program code by reading the executable program code stored in the memory, so as to implement the gene regulatory network inference method of spatiotemporal transcriptome data.

[0051] In another aspect, the present disclosure provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the gene regulatory network inference method for spatiotemporal transcriptome data.

[0052] Additional aspects and advantages of the present disclosure will be given in part in the description that follows and, in part, will be obvious from the description that follows, or will be learned through practice of the present disclosure. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] The above and / or additional aspects and advantages of the present disclosure will become apparent and readily understood from the description of the embodiments in conjunction with the following drawings, in which:

[0054] FIG1 shows a schematic flow chart of a method for inferring gene regulatory networks from spatiotemporal transcriptome data according to an embodiment of the present disclosure;

[0055] Figure 2 shows a framework diagram of a gene regulatory network inference method for spatiotemporal transcriptome data;

[0056] FIG3 shows the preliminary results of cell dimensionality reduction clustering in an embodiment of the present disclosure.

[0057] Detailed Description of the Invention

[0058] The embodiments of the present disclosure are described in detail below. The embodiments described below are exemplary and are only used to explain the present disclosure, and should not be understood as limiting the present disclosure.

[0059] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features being referred to. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one such feature. Throughout the present disclosure, "plurality" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0060] The endpoints of the ranges and any values ​​disclosed herein are not limited to the precise ranges or values, and these ranges or values ​​should be understood to include values ​​close to these ranges or values. For numerical ranges, the endpoints of each range, the endpoints of each range and individual point values, and the individual point values ​​can be combined with each other to obtain one or more new numerical ranges, which should be considered to be specifically disclosed herein.

[0061] In order to make the present disclosure more easily understood, certain technical and scientific terms are specifically defined below. Unless otherwise clearly defined elsewhere in this document, all other technical and scientific terms used herein have the meaning commonly understood by those skilled in the art to which the present disclosure belongs.

[0062] In this document, the terms “include” or “comprising” are open expressions, that is, including the contents specified in the present disclosure, but not excluding other contents.

[0063] As used herein, the terms "optionally," "optional," or "optionally" generally mean that the subsequently described event or circumstance may but need not occur, and that the description includes instances where the event or circumstance occurs and instances where it does not.

[0064] In the present disclosure, the term "gene expression matrix" refers to the data results obtained using spatial transcriptomics technology. In the matrix, rows represent the expression of a gene under different environmental conditions or at different time points, and columns represent the expression of all genes under different conditions or samples (such as tissues, experimental conditions, treatment factors, etc.). The data in each grid represents the expression level of a specific gene in a specific sample.

[0065] In this disclosure, the term "structural equation modeling" refers to a method for constructing, estimating, and testing causal models. Models include both observable explicit variables and potentially unobservable latent variables. Structural equation modeling can replace methods such as multiple regression, path analysis, factor analysis, and covariance analysis to clearly analyze the impact of individual indicators on the overall population and the relationships between individual indicators.

[0066] In this disclosure, the term "variational autoencoder" refers to an autoencoder, which is a model that discovers some latent state (incomplete, sparse, denoised, contracted) representation of the data. More specifically, the input data is converted into an encoded vector where each dimension represents an attribute learned from the data. The most important thing is that the encoder outputs a single value for each encoded dimension, and the decoder then takes these values ​​and tries to recreate the original input. The variational autoencoder (VAE) provides a probabilistic way to describe observations in the latent space. Therefore, instead of building an encoder that outputs a single value to describe each hidden state attribute, the encoder describes the probability distribution of each hidden attribute.

[0067] In this disclosure, the term "dimensionality reduction" refers to the attempt to compress the dimensions while preserving the distribution information as much as possible. We can think of it as data compression or feature selection.

[0068] In this disclosure, the term "adjacency matrix" is logically structured into two parts: V and E sets, where V represents vertices and E represents edges. Therefore, a one-dimensional array stores all vertex data in a graph, while a two-dimensional array stores data on the relationships between vertices (edges or arcs). This two-dimensional array is called the adjacency matrix. Adjacency matrices are further divided into directed graph adjacency matrices and undirected graph adjacency matrices.

[0069] In this disclosure, the term "Gaussian distribution" is also called normal distribution.

[0070] In this disclosure, the term "normalization" refers to scaling the feature value to the interval [0, 1] using the maximum and minimum values ​​of a feature (which can be understood as a data set), and scaling the features of each column using the min-max function.

[0071] In this disclosure, the term "characteristic matrix" generally refers to the eigenvalue of a matrix. For example, if A is an n-order square matrix, if there exists a number m and a non-zero n-dimensional column vector x such that Ax=mx, then m is called an eigenvalue or eigenvalue of the matrix A.

[0072] In this disclosure, the term "hyperparameter" refers to parameters that determine the model. Different models have different hyperparameters, and hyperparameters are generally variables determined empirically. In deep learning, hyperparameters include: learning rate, number of iterations, number of layers, number of neurons in each layer, etc. In this disclosure, the learning rate refers to the step size that determines each iteration of the gradient descent algorithm.

[0073] In this disclosure, the term "loss function" is an operational function used to measure the degree of difference between the model's predicted value f(x) and the true value Y. It is a non-negative real-valued function, usually represented by L(Y, f(x)). The smaller the loss function, the better the robustness of the model. The loss function is mainly used in the training phase of the model. After each batch of training data is fed into the model, the predicted value is output through forward propagation, and then the loss function calculates the difference between the predicted value and the true value, which is the loss value. After obtaining the loss value, the model updates each parameter through back propagation to reduce the loss between the true value and the predicted value, so that the predicted value generated by the model moves closer to the true value, thereby achieving the purpose of learning.

[0074] In this disclosure, the term "structural equation modeling" refers to a method for building, estimating, and testing causal models. It can replace methods such as multiple regression, path analysis, factor analysis, and covariance analysis to clearly analyze the impact of individual indicators on the overall population and the relationships between individual indicators.

[0075] The element aij in the gene expression matrix represents the expression level of gene j in cell i. The element bij in the cell position matrix is ​​calculated by calculating the normalized distance between cells i and j using their coordinates. This is done by calculating the normalized distance between the two cells as 1-ij. In other words, cells with greater distances have smaller bij values, while cells with closer distances have larger bij values. Adjacency is determined by the magnitude of bij. Generally, cells i and j are considered adjacent when bij > 0.5.

[0076] A gene regulatory network inference method for spatiotemporal transcriptome data

[0077] Starting from graph neural networks, this paper uses gene expression matrices and cell location information from spatiotemporal data to infer gene regulatory networks. This paper designs a graph variational network based on graph neural networks, using the adjacency matrix as a trainable parameter and reconstructing input features through a variational autoencoder to identify regulatory relationships between genes.

[0078] According to some specific embodiments of the present disclosure, the present disclosure provides a method for inferring gene regulatory networks from spatiotemporal transcriptome data, as shown in FIG1 , the method comprising:

[0079] S100, using the spatiotemporal transcriptome platform to obtain the gene expression matrix and cell location matrix of the sample;

[0080] S110, using cells as feature dimensions and genes and cells as node dimensions, merging the gene expression matrix and the cell position matrix to obtain a node feature matrix X with genes and cells as nodes and cells as features;

[0081] S120, initializing an adjacency matrix A representing the connection information of the graph nodes;

[0082] S130, inputting the initialized adjacency matrix A and node feature matrix X into the encoder part of the autoencoder, performing dimensionality reduction processing by the encoder, and extracting the embedding vector;

[0083] S140, inputting the embedding vector into the decoder part of the autoencoder to obtain a reconstructed node feature matrix Y;

[0084] S150 . Obtain a trained adjacency matrix A1 based on the node feature matrix X and the reconstructed node feature matrix Y.

[0085] According to some specific embodiments of the present disclosure, the present disclosure provides a method for inferring gene regulatory networks from spatiotemporal transcriptome data, the method comprising:

[0086] (1) Using the spatiotemporal transcriptome platform to obtain the gene expression matrix and cell location matrix of the sample, cells are used as the feature dimension, genes and cells are used as the node dimension, and the two matrices are merged to obtain the node feature matrix X with genes + cells as nodes and cells as features;

[0087] (2) Randomly initialize the adjacency matrix A that represents the connection information of the graph nodes;

[0088] (3) The constructed node adjacency matrix A and node feature matrix X are input into the encoder part of the graph variational autoencoder. When input into the encoder, the hyperparameter n is introduced. channel The node feature matrix X of the two-dimensional tensor is expanded into a three-dimensional tensor, and finally the encoder is used for dimensionality reduction to extract the embedding vectors of genes and cells respectively;

[0089] (4) The embedded vector is input into the decoder of the graph variational autoencoder. During this period, the feature matrix of the node is reconstructed by optimizing the loss function, introducing the re-parameterization and Gumbel Softmax techniques;

[0090] (5) Extract the adjacency matrix A obtained after training to obtain the gene regulatory network.

[0091] According to the specific implementation scheme of the present disclosure, the main framework of this algorithm (Figure 2) is as follows: the gene expression matrix and the cell position matrix are obtained from the spatiotemporal transcriptome, and the above two matrices are merged into a node feature matrix X (taking cells as the feature dimension, genes and cells as the node dimension, and merging the two matrices vertically to obtain a feature matrix with genes + cells as nodes and cells as features) as the input of the graph variational autoencoder. The adjacency matrix A representing the connection information of the graph nodes is randomly initialized, which contains three parts: gene-gene adjacency relationship (i.e., the gene regulatory network to be inferred), cell-cell adjacency relationship, and gene-cell adjacency relationship. The constructed node adjacency matrix and feature matrix are input into the encoder part of the graph variational autoencoder. After the encoder dimensionality reduction processing, the embedding vectors of genes and cells are extracted respectively, and the embedding vectors are input into the decoder part of the graph variational autoencoder to reconstruct the feature matrix of the node. During the training process, the adjacency matrix will aggregate the neighboring node information of each node in the graph to the node. When the graph variational autoencoder reconstructs the node feature matrix, the adjacency matrix A as a trainable parameter will find the optimal aggregation relationship of the node features, thereby inferring the gene regulatory network.

[0092] According to specific embodiments of the present disclosure, the adjacency matrix A comprises three components: gene-gene adjacency, cell-cell adjacency, and gene-cell adjacency. All three components can be randomly initialized, especially the gene-gene adjacency. The cell-cell and gene-cell adjacency components can also be determined via input. These two components are also updated as parameters during training. While determining these adjacency components via input does not affect the final result, it may affect the speed of parameter updates.

[0093] According to a specific embodiment of the present disclosure, the spatiotemporal transcriptome platform includes at least one selected from the following: the Biotron S1000 spatial transcriptome sequencing platform, the STOmics spatial transcriptome sequencing platform, and the 10xgenomics spatial transcriptome sequencing platform, and the samples include tissue samples selected from animals and plants.

[0094] According to a specific embodiment of the present disclosure, both the encoder part and the decoder part of the graph variational autoencoder use graph convolutional layers for feature aggregation and reconstruction.

[0095] According to a specific embodiment of the present disclosure, the gene regulatory network adopts 3 hidden layers and is inferred based on a neural network model.

[0096] According to a specific embodiment of the present disclosure, after the adjacency matrix A is normalized, it and the node feature matrix X jointly determine the Gaussian distribution Z obeyed by the embedded features.

[0097] According to a specific embodiment of the present disclosure, the optimization of the loss function includes introducing an L1 regularization term for the adjacency matrix, and the optimization range includes the loss of feature reconstruction, the difference between the embedded feature distribution and the Gaussian distribution, and the L1 regularization of the adjacency matrix.

[0098] According to a more specific embodiment of the present disclosure, a gene regulatory network inference method for spatiotemporal transcriptome data is provided. Both the encoder and decoder parts use classic graph convolutional layers for feature aggregation and reconstruction. The method includes:

[0099] (1) Encoder part

[0100] For the input node feature X and adjacency matrix A, the node feature aggregation is performed through the graph convolutional network (GCN) to obtain the expected feature distribution and variance The Gaussian distribution Z that the embedded features obey is determined by μ and logσ, where It is a normalized adjacency matrix. The present disclosure uses a directed acyclic graph to represent the gene regulatory network. In the present disclosure, the node feature X input to the encoder is a two-dimensional tensor, with rows representing each node and columns representing high-dimensional features of each node. In gene regulatory network inference, X is a gene expression matrix, with rows representing genes and columns representing cells. Assuming that the gene expression matrix is ​​a two-dimensional tensor consisting of n genes and m cells, that is, X = X(n,m). In order to make the encoder utilize the node feature X to the greatest extent, we expand the two-dimensional tensor X into a three-dimensional tensor, that is, X = X(n,m,1). When X passes through the encoder (in the present disclosure, a graph convolution layer is used for node feature aggregation and dimensionality transformation), the third dimension is dimensionally increased or decreased, and the hyperparameter n is introduced. channel , that is, when X passes through the encoder, the dimension of X becomes X=X(n,m,n channel ). When passing through the last layer of the encoder, the dimension of X becomes X = X(n,m,1), and the third dimension is compressed, and the node features can be converted back to a two-dimensional tensor to calculate the Gaussian distribution that the features need to obey. In this disclosure, n is selected. channel =128.

[0101] (2) Decoder

[0102] The present invention reconstructs node features, while the traditional graph variational autoencoder reconstructs the adjacency matrix A of the graph, using the embedding vector Z multiplied by its own transpose and then passing through the Sigmoid activation function, that is, in It is a reconstructed graph adjacency matrix, so the decoder part is different from the traditional graph variational autoencoder. The embedding vector Z obtained by the encoder is passed through the graph convolution layer (i.e., the decoder) to obtain a reconstructed feature matrix consistent with the input feature dimension. In the decoder part, as in the encoder part, the dimension of the node feature tensor X is converted to three dimensions for decoding in order to maximize the use of node features. Because the encoder extracts features from genes + cells as graph nodes, the obtained embedding vector contains both gene embeddings and cell embeddings, that is, the embedding vector obtained by splicing the cell embeddings and the gene embeddings. By keeping the node order of this embedding vector consistent with the feature matrix input to the encoder, those skilled in the art can extract cell embeddings from the embeddings for preliminary cell dimensionality reduction and clustering analysis.

[0103] Since the inventors of this application use a graph convolution layer as the network architecture, that is, node feature extraction and aggregation are performed through the adjacency matrix A, after setting A as a trainable parameter, the encoder-decoder architecture will learn the optimal adjacency matrix A by reconstructing the node features after aggregation by the adjacency matrix A, thereby better reflecting the connection relationship between nodes. In this patent, the gene-gene relationship contained in the critical matrix, that is, the regulatory relationship between transcription factors and target genes, is characterized. By extracting the relationship pairs revealed in the adjacency matrix A, the gene regulatory relationship is obtained.

[0104] (3) Loss function

[0105] In this disclosure, the inventors introduce the L1 regularization term of the adjacency matrix based on the loss function of the classic graph variational autoencoder, and the loss function is L=E q(Z|X,A) [log p(A|Z)] - KL[q(Z|X,A)||p(z)] + αL1(A), where KL is the Kullback-Leibler divergence and α is a tunable hyperparameter. This loss function represents the loss of feature reconstruction, the difference between the embedded feature distribution and the Gaussian distribution, and L1 regularization of the adjacency matrix.

[0106] By simultaneously optimizing the three loss functions during training, we can ensure that (I) the reconstructed node feature tensor is as consistent as possible with the input feature tensor distribution; (II) the embedded feature distribution extracted by the variational autoencoder is as close to a Gaussian distribution as possible, where each input feature corresponds to a Gaussian distribution, so that the decoder can sample from this Gaussian distribution; and (III) the adjacency matrix representing the gene regulatory network exhibits a certain sparsity, which can be regularized by adjusting the hyperparameter α.

[0107] The new neural network parameters can be updated through the loss function composed of the above three parts, namely, reconstructing features, embedding vectors to obey Gaussian distribution, and normalizing the adjacency matrix.

[0108] (4) Heavy parameters and Gumbel Softmax technique

[0109] After the encoder obtains the expected μ and variance logσ of the feature distribution, the expected and variance determine the Gaussian distribution Z that the feature follows. The output Z is a random node. Reparameterization transforms the random Z into a deterministic node while introducing another node ∈ to create randomness. Therefore, after obtaining the expected and variance, the Gaussian distribution that the feature follows is Z = μ + ∈·logσ. This approach removes the sampling process from the computational graph, allowing the entire graph to calculate gradient updates. ∈ is not updated during the computational process and is treated as a parameter with no weight change.

[0110] For the n-dimensional probability vector π, the corresponding discrete random variable x π After adding Gumbel noise and sampling, we get: x π =argmax(log(π i )+G i )

[0111] Among them, G i is an independent and identically distributed standard Gumbel distributed random variable, and the probability density function of the standard Gumbel distribution is: Because argmax is not differentiable, softmax is used instead of argmax, which is the Gumbel Softmax technique. i It can be obtained from the uniform distribution by inverting the Gumbel distribution, that is, G i =-log(-log(U i )), where U i ~U(0,1) is a uniform distribution.

[0112] Introducing the Gumbel softmax technique adds randomness to simulate distributed sample generation, increasing the exploration space and model robustness. Gumbel softmax operates on the embedding vectors generated by the encoder, effectively approximating the probability distribution of the network output. It also performs sampling and derivatives, making the Gaussian-distributed embeddings more robust before passing through the decoder.

[0113] The present disclosure uses intermediate variables to perform cell dimensionality reduction and clustering while embedding information, which can be preliminarily used for cell dimensionality reduction and clustering. The present disclosure performs gene regulatory network inference on spatiotemporal transcriptome data, especially combines cell position information to perform gene regulatory network inference, makes full use of spatiotemporal information, minimizes information loss in the gene regulatory network inference process, and thus improves the accuracy of inference.

[0114] Example 1 Evaluation test of algorithm model

[0115] Mouse embryos with a developmental period of E9.5 were selected, and spatial transcriptome sequencing was performed using BGI Stereo-seq to obtain an E9.5 dataset. The dataset was read using the Python third-party library Scanpy to obtain the expression matrix and cell position matrix. The data were tested using different algorithm models, and the test results were compared with the regulatory relationships obtained in the experiment to score the training results. The scoring standard used was EPR, and a higher score indicated that the training result was closer to the experimental result.

[0116] The experimental results are shown in Table 1, which show that the model of the present invention outperforms SCENIC and DeepSEM. The SCENIC model extracts each gene, takes each target gene as output, and takes the transcription factor as input, constructs n random forest trees (n is the number of genes), and calculates the importance score between each transcription factor and the gene, and finally obtains the transcription factor-gene co-expression module, and finally deletes the gene relationship below the threshold; the DeepSEM model refers to embedding the structural equation model into the variational autoencoder to achieve the prediction of the interaction between variables; Ours, that is, the model of the present invention, refers to using the graph variational autoencoder as a framework, by inputting the gene expression matrix and the cell position matrix, and training the adjacency matrix A that represents the gene regulatory relationship in the process of reconstructing the input features.

[0117] Table 1

[0118] The specific experimental process is as follows:

[0119] 1) Obtain spatial transcriptome data of mouse embryos at developmental stage E9.5;

[0120] 2) Input the dataset obtained in step 1) into Scanpy to obtain the expression matrix and cell position matrix of the spatial transcriptome;

[0121] 3) Merge the expression matrix and the cell position matrix vertically to obtain the node feature matrix X with genes + cells as nodes and cells as features;

[0122] 4) Randomly initialize the adjacency matrix A that represents the connection information of the graph nodes;

[0123] 5) The constructed adjacency matrix and node feature matrix are input into the encoder part of the graph variational autoencoder. After the encoder performs dimensionality reduction processing, the embedding vectors of genes and cells are extracted respectively. The embedding vectors are input into the decoder part of the graph variational autoencoder to reconstruct the node feature matrix Y;

[0124] 6) During the training process, the adjacency matrix A aggregates the neighboring node information of each node in the graph to that node. When the graph variational autoencoder reconstructs the node feature matrix X, the adjacency matrix A, as a trainable parameter, searches for the optimal aggregation relationship of the node features, thereby training the adjacency matrix A and obtaining the trained adjacency matrix A1.

[0125] 7) Extract gene regulatory relationships based on trained A1;

[0126] 8) Compare the extracted gene regulatory relationships with the regulatory relationships obtained experimentally to obtain a training score.

[0127] At the same time, the intermediate generated feature distribution is subjected to PCA dimensionality reduction and then visualized to preliminarily obtain cell dimensionality reduction clustering information, as shown in Figure 3, which shows that the cell dimensionality reduction clustering results obtained using the disclosed model are accurate.

[0128] Based on spatiotemporal transcriptome data, this paper proposes and designs a gene regulatory network inference method based on graph neural network, which provides a reliable solution for inferring gene regulatory networks for spatiotemporal transcriptome data.

[0129] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present disclosure. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine different embodiments or examples described in this specification and features of different embodiments or examples, unless they are mutually inconsistent.

[0130] Although the embodiments of the present disclosure have been shown and described above, it is understood that the above embodiments are illustrative and are not to be construed as limitations on the present disclosure. A person skilled in the art may change, modify, replace and vary the above embodiments within the scope of the present disclosure.

Claims

1. A gene regulatory network inference method based on spatiotemporal transcriptome data, characterized in that: The method comprises: (A) Gene expression matrix and cell location matrix of samples were obtained using spatiotemporal transcriptome platform; (B) Taking cells as feature dimensions and genes and cells as node dimensions, the gene expression matrix and the cell position matrix are merged to obtain a node feature matrix X with genes and cells as nodes and cells as features; (C) Initialize the adjacency matrix A that represents the connection information of the graph nodes; (D) inputting the initialized adjacency matrix A and node feature matrix X into the encoder part of the autoencoder, and extracting the embedding vector after the encoder performs dimensionality reduction processing; (E) inputting the embedding vector into the decoder part of the autoencoder to obtain a reconstructed node feature matrix Y; (F) Based on the node feature matrix X and the reconstructed node feature matrix Y, obtain the trained adjacency matrix A1.

2. The method according to claim 1, characterized in that In step (B), the direction of matrix merging is vertical merging.

3. The method according to claim 1 or 2, characterized in that: The sample includes a tissue sample selected from animals and plants.

4. The method according to any one of claims 1 to 3, characterized in that: In step (C), the adjacency matrix A includes three parts: gene-gene adjacency relationship, cell-cell adjacency relationship, and gene-cell adjacency relationship.

5. The method according to any one of claims 1 to 4, characterized in that: In step (D), the autoencoder includes a graph variational autoencoder.

6. The method according to claim 5, characterized in that The encoder part of the graph variational autoencoder uses a graph convolutional layer for feature aggregation, and the decoder part of the graph variational autoencoder uses a graph convolutional layer for reconstruction.

7. The method according to any one of claims 1 to 6, characterized in that: In step (D), when the adjacency matrix A and the node feature matrix X are input into the encoder part, the hyperparameter n is introduced channel The node feature matrix X of the two-dimensional tensor is expanded to a three-dimensional tensor, and after the dimensionality reduction processing by the encoder, the embedding vectors of genes and cells are extracted respectively.

8. The method according to any one of claims 1 to 7, characterized in that: In step (E), the reconstructed node feature matrix Y is obtained by the following method: The embedding vector is input into the decoder part of the graph variational autoencoder, and the feature matrix Y of the node is reconstructed by optimizing the loss function, introducing the heavy parameters and Gumbel Softmax technique.

9. The method according to claim 8, characterized in that The optimization of the loss function includes introducing an L1 regularization term for the adjacency matrix. Furthermore, the optimization of the loss function includes optimization in the process of reconstructing the feature matrix Y and optimization in the process of making the embedded vector obey a Gaussian distribution.

10. The method according to any one of claims 1 to 9, characterized in that: After step (F), the method further includes step (G), obtaining gene regulatory relationships based on the trained adjacency matrix A1.

11. A system for inferring gene regulatory networks based on spatiotemporal transcriptome data, characterized in that: include: A data acquisition module, which is used to acquire a gene expression matrix and a cell location matrix of a sample using a spatiotemporal transcriptome platform; A matrix merging module, the matrix merging module is connected to the data acquisition module, the matrix merging module takes cells as feature dimensions, genes and cells as node dimensions, and vertically merges the gene expression matrix and the cell position matrix to obtain a node feature matrix X with genes and cells as nodes and cells as features; An initialization module, the initialization module is connected to the matrix merging module, and the initialization module is used to initialize the adjacency matrix A representing the connection information of the graph nodes; A dimensionality reduction processing module, which is connected to the initialization module, and is used to input the initialized adjacency matrix A and node feature matrix X into the encoder part of the autoencoder, and extract the embedding vector after the encoder dimensionality reduction processing; A reconstruction module, the reconstruction module is connected to the dimension reduction processing module, and the reconstruction module is used to input the embedding vector into the decoder part of the autoencoder to obtain a reconstructed node feature matrix Y; A feature training module, the feature training module is connected to the reconstruction module, and the feature training module is used to obtain a trained adjacency matrix A1 based on the node feature matrix X and the reconstructed node feature matrix Y.

12. The system according to claim 11, wherein: In the matrix merging module, the direction of matrix merging is vertical merging.

13. The system according to claim 11 or 12, wherein: The sample includes a tissue sample selected from animals and plants.

14. The system according to any one of claims 11 to 13, characterized in that: The adjacency matrix A includes three parts: gene-gene adjacency relationship, cell-cell adjacency relationship, and gene-cell adjacency relationship.

15. The system according to any one of claims 11 to 14, characterized in that: The autoencoder comprises a graph variational autoencoder.

16. The system according to claim 15, characterized in that The encoder part of the graph variational autoencoder uses a graph convolutional layer for feature aggregation, and the decoder part of the graph variational autoencoder uses a graph convolutional layer for reconstruction.

17. The system according to any one of claims 11 to 16, characterized in that: In the dimensionality reduction processing module, when the adjacency matrix A and the node feature matrix X are input into the encoder part, the hyperparameter n is introduced channel The node feature matrix X of the two-dimensional tensor is expanded to a three-dimensional tensor, and after the dimensionality reduction processing by the encoder, the embedding vectors of genes and cells are extracted respectively.

18. The system according to any one of claims 11 to 17, characterized in that: In the reconstruction module, the reconstructed node feature matrix Y is obtained by the following method: The embedding vector is input into the decoder part of the graph variational autoencoder, and the feature matrix Y of the node is reconstructed by optimizing the loss function, introducing the heavy parameters and Gumbel Softmax technique.

19. The system according to claim 18, characterized in that The optimization of the loss function includes introducing an L1 regularization term for the adjacency matrix. Furthermore, the optimization of the loss function includes optimization in the process of reconstructing the feature matrix Y and optimization in the process of making the embedded vector obey a Gaussian distribution.

20. The system according to any one of claims 11 to 19, characterized in that: After the feature training module, a data extraction module is further included, and the data extraction module is connected to the feature training module. The data extraction module obtains the gene regulation relationship based on the trained adjacency matrix A1.

21. An electronic device, characterized in that: Including memory and processor; The processor runs a program corresponding to the executable program code by reading the executable program code stored in the memory, so as to implement the gene regulatory network inference method of spatiotemporal transcriptome data described in any one of claims 1 to 10.

22. A computer-readable storage medium storing a computer program, characterized in that: When the program is executed by a processor, the gene regulatory network inference method of spatiotemporal transcriptome data according to any one of claims 1 to 10 is implemented.

Images