A dimension reduction method and system for spatial transcriptome data

Abstract

The application discloses a dimension reduction method and system for spatial transcriptome data, belongs to the technical field of data analysis, and can solve the problem that the existing dimension reduction method has high calculation complexity and cannot simultaneously perform spatial dimension reduction on multiple slices. The method comprises the following steps: step 1, acquiring gene expression matrices and spatial coordinate matrices of multiple slices in spatial transcriptome data; step 2, determining initial dimension reduction parameters and spatial dependencies between spots in each slice according to the gene expression matrices and the spatial coordinate matrices; step 3, constructing a dimension reduction model according to the gene expression matrices, the initial dimension reduction parameters and the spatial dependencies between the spots in the slices; and step 4, estimating parameters of the dimension reduction model by using an iterative optimization algorithm to obtain a final low-dimensional representation of each slice. The application is used for dimension reduction of spatial transcriptome data.

Description

A method and system for dimensionality reduction of spatial transcriptome data Technical Field

[0001] This invention relates to a method and system for dimensionality reduction of spatial transcriptome data, belonging to the field of data analysis technology. Background Technology

[0002] Existing literature contains numerous dimensionality reduction methods, primarily categorized into parametric and nonparametric methods. Parametric methods include Principal Component Analysis (PCA) and Zero-Inflation Factor Analysis (ZIFA), while nonparametric methods include T-Stochastic Neighbor Embedding (tSNE), Uniform Manifold Approximation and Projection (UMAP), and Single-Cell Variational Inference (scVI). PCA is the most common dimensionality reduction method, widely used in many software workflows; tSNE and UMAP are two widely used nonparametric dimensionality reduction methods for visualization; and scVI is a deep generative representation learning method for single-cell RNA sequencing data analysis. However, none of these methods consider the spatial characteristics of Spatially Resolved Transcriptoms (SRT) data when generating low-dimensional embeddings.

[0003] In SRT data analysis, gene expression at adjacent spatial locations exhibits "similarity" due to shared cellular microenvironments. Recently proposed SpatialPCA and non-negative spatial factorization (NSF) utilize Gaussian kernels and location-based Gaussian process priors, respectively, to achieve spatially perceptual dimensionality reduction. To reduce computational burden, SpatialPCA employs a low-rank approximation strategy, achieving O(n^2) computational efficiency. 2 The computational complexity is calculated as n, where n is the number of spatial locations. NSF employs a sparse Gaussian process, whose computational complexity is proportional to the cube of the number of induction points, but it is not suitable for spatial locations from multiple slices. With the improvement of spatial resolution in SRT technology, the number of spatial locations increases significantly, requiring multiple slices to generate spatial maps or reconstruct spatiotemporal transcriptomic maps of the entire organ. However, existing dimensionality reduction methods have high computational complexity and cannot simultaneously perform spatial dimensionality reduction on multiple slices. Summary of the Invention

[0004] This invention provides a method and system for dimensionality reduction of spatial transcriptome data, which can solve the problems of high computational complexity and inability to perform spatial dimensionality reduction on multiple slices simultaneously by existing dimensionality reduction methods.

[0005] On one hand, the present invention provides a method for dimensionality reduction of spatial transcriptome data, the method comprising:

[0006] Step 1: Obtain the gene expression matrix and spatial coordinate matrix of multiple slices from the spatial transcriptome data;

[0007] Step 2: Determine the initial dimensionality reduction parameters and the spatial dependence between spots in each slice based on the gene expression matrix and the spatial coordinate matrix;

[0008] Step 3: Construct a dimensionality reduction model based on the gene expression matrix, the initial dimensionality reduction parameters, and the spatial dependence between spots in the slice;

[0009] Step 4: Use an iterative optimization algorithm to estimate the parameters of the dimensionality reduction model and obtain the final low-dimensional representation of each slice.

[0010] Optionally, step 2 specifically includes:

[0011] Step 21: Based on the gene expression matrix, use principal component analysis to determine the initial dimensionality reduction parameters and the initial low-dimensional representation of each slice;

[0012] Step 22: Determine the spatial dependency between spots in each slice based on the spatial coordinate matrix and the initial low-dimensional representation of each slice.

[0013] Optionally, the spatial dependency between spots in the slice is specifically as follows:

[0014] Given low-dimensional representations of all spot samples other than the target spot in the slice, the low-dimensional representation of the target spot follows a normal distribution; wherein the mean of the normal distribution is the average of the low-dimensional representations of the neighbors of the target spot.

[0015] Optionally, step 1 specifically includes:

[0016] Step 11: Obtain the hypervariable genes and spatial coordinate matrix of each slice in the spatial transcriptome data to be analyzed, and integrate the hypervariable genes of each slice to obtain the integrated hypervariable genes.

[0017] Step 12: Obtain the gene expression matrix of each slice based on the integrated hypervariable genes.

[0018] Optionally, step 12 specifically includes:

[0019] Step 121: Extract the counting matrix of the integrated hypervariable genes in each slice;

[0020] Step 122: Regularize the counting matrix of each slice to obtain the gene expression matrix of each slice.

[0021] Optionally, the dimensionality reduction model is:

[0022] The gene expression matrix of each slice is the sum of three terms: the sample-specific intercept term of the slice, the product of the sample shared load matrix and the low-dimensional representation of the slice, and the random noise of the slice.

[0023] Optionally, the dimensionality reduction model is specifically as follows:

[0024] x mi =α m +Wv mi +∈ mi ;

[0025] Where, x mi ∈R p Let α be the gene expression vector of the i-th spot in the m-th slice. m ∈R p For the sample-specific intercept term of the m-th slice, W∈R p×q For the sample shared load matrix, v mi ∈R q For the low-dimensional representation of the i-th spot in the m-th slice, ∈ mi ~N(0,Λ m ) represents the normally distributed random noise of the i-th spot in the m-th slice, where Λ m =diag(λ mi ,…,λ mp ) is a diagonal matrix.

[0026] Optionally, the dimensionality reduction model is:

[0027] Given the Poisson rate of the target spot in a given slice, the gene expression of the target spot follows a Poisson distribution;

[0028] Furthermore, the logarithm of the Poisson rate of the target spot is the sum of three terms: the sample-specific intercept term of the slice, the product of the sample shared load vector and the low-dimensional representation of the slice, and the random noise of the slice.

[0029] Optionally, the dimensionality reduction model is specifically as follows:

[0030]

[0031]

[0032] in, Let a be the expression level of gene j in the i-th spot of slice m. mi w is the regularization factor for the i-th spot in slice m. j ∈R q Let α be the sample shared payload vector corresponding to gene j. mj f is the sample-specific intercept term corresponding to gene j in slice m. mij v is the Poisson rate of gene j in the i-th spot of slice m; mi ∈R q Let m be the low-dimensional representation of the i-th spot in slice m, ∈ mij ~N(0,λ) mj ) represents the random noise of the normally distributed gene j in the i-th spot of slice m.

[0033] Optionally, the iterative optimization algorithm is the ICM-EM algorithm.

[0034] Optionally, the regularization process is logarithmic regularization or negative binomial regression regularization.

[0035] On the other hand, the present invention provides a dimensionality reduction system for spatial transcriptome data, the system comprising:

[0036] The data acquisition module is used to acquire the gene expression matrix and spatial coordinate matrix of multiple slices in the spatial transcriptome data;

[0037] The parameter determination module is used to determine the initial dimensionality reduction parameters and the spatial dependence between spots in each slice based on the gene expression matrix and the spatial coordinate matrix.

[0038] The model building module is used to construct a dimensionality reduction model based on the gene expression matrix, the initial dimensionality reduction parameters, and the spatial dependence between spots in the slice.

[0039] The iterative optimization module is used to estimate the parameters of the dimensionality reduction model using an iterative optimization algorithm, and obtain the final low-dimensional representation of each slice.

[0040] The beneficial effects that this invention can produce include:

[0041] The spatial transcriptome data dimensionality reduction method provided by this invention, denoted as ProFAST, can effectively estimate low-dimensional representations across multiple slices, capturing inherent biological effects while considering local expression similarities caused by shared cellular microenvironments. Compared to existing technologies, ProFAST allows for simultaneous spatial dimensionality reduction across multiple slices while modeling counting properties. Furthermore, ProFAST uses conditional autoregression to model local spatial dependencies, achieving O(n) computational complexity, making it suitable for analyzing multi-slice high-resolution SRT data. Attached Figure Description

[0042] Figure 1 is a flowchart of the spatial transcriptome data dimensionality reduction method provided in an embodiment of the present invention;

[0043] Figure 2 is a schematic diagram of the dimensionality reduction principle of spatial transcriptome data provided in the embodiments of the present invention. Detailed Implementation

[0044] The present invention will be described in detail below with reference to the embodiments, but the present invention is not limited to these embodiments.

[0045] This invention provides a dimensionality reduction method for spatial transcriptome data, denoted as ProFAST; as shown in Figures 1 and 2, the method includes:

[0046] Step 1: Obtain the gene expression matrix and spatial coordinate matrix of multiple slices in the spatial transcriptome data.

[0047] Step 1 specifically includes:

[0048] Step 11: Obtain the hypervariable genes and spatial coordinate matrix of each slice in the spatial transcriptome data to be analyzed, and integrate the hypervariable genes of each slice to obtain the integrated hypervariable genes.

[0049] Step 12: Obtain the gene expression matrix for each slice based on the integrated hypervariable genes.

[0050] Step 12 specifically includes:

[0051] Step 121: Extract the counting matrix of the integrated hypervariable genes in each slice;

[0052] Step 122: Regularize the counting matrix of each slice to obtain the gene expression matrix of each slice.

[0053] In this embodiment of the invention, assuming the number of slices is M, the gene expression matrix can be denoted as X. m The spatial coordinate matrix can be denoted as P. m Gene expression matrix X, m = 1, ..., M. m It is an n m A matrix of size ×p, where nm Let p be the sample size of the m-th slice, and p be the number of hypervariable genes. The gene expression matrix can be a UMI count matrix of gene expression or a regularized expression matrix. Regularization methods can include logarithmic regularization or negative binomial regression regularization. Spatial coordinate matrix P m It is an n m A 2×2 matrix, where the i-th row represents the spatial coordinates of the i-th spot in the m-th slice. In the context of spatial transcriptomics, a spot refers to a discrete location in a tissue sample where gene expression data is collected. These spots are typically identified by physical coordinates or imaging techniques such as microscopy. From a practical perspective, a spot can be considered a cellular unit, representing a single cell or a group of cells in a tissue.

[0054] For example, the data to be analyzed may include UMI count matrices and spatial coordinate information (i.e., spatial coordinate matrix P) of multiple spatial transcriptome slices obtained from genomics databases (such as GEO). m The *FindSVGs* function in the R package DR.SC is used to select hypervariable genes for each slice, or the *FindVariableFeatures* function in the R package Seurat is used to select hypervariable genes for each slice. The hypervariable genes from each slice are then integrated according to the sorting rules to obtain the integrated hypervariable gene matrix. The UMI counts on the integrated hypervariable genes are then regularized using the *NormalizeData* function or the *SCTransform* function in the R package Seurat to obtain the regularized gene expression matrix X. m .

[0055] Step 2: Determine the initial dimensionality reduction parameters and the spatial dependence between spots in each slice based on the gene expression matrix and spatial coordinate matrix;

[0056] Step 2 specifically includes:

[0057] Step 21: Based on the gene expression matrix X m Principal component analysis (PCA) was used to determine the initial dimensionality reduction parameters and the initial low-dimensional representation of each slice.

[0058] Step 22: Determine the spatial dependencies between spots in each slice based on the spatial coordinate matrix and the initial low-dimensional representation of each slice.

[0059] The spatial dependence between spots in a slice is specifically as follows:

[0060] Given the low-dimensional representations of all spot samples outside the target spot in a given slice, the low-dimensional representation of the target spot follows a normal distribution; where the mean of the normal distribution is the average of the low-dimensional representations of the target spot's neighbors.

[0061] Specifically, based on the gene expression matrix X m and spatial coordinate matrix P m The initial dimensionality reduction parameters (i.e., W, Ψ) were determined using principal component analysis (PCA). m ,Λ m ) and the initial low-dimensional representation of each slice Where v mi It is a low-dimensional representation of the i-th spot in the m-th slice.

[0062] To capture the spatial dependencies between spots in each slice, v mi Modeling is performed using continuous multivariate hidden Markov random fields. Specifically, assume v mi It comes from a CAR (Conditional Autoregressive) model.

[0063] The spatial dependence between spots in the slice is shown in equation (1):

[0064]

[0065] Where the subscript [n m ]\i represents all samples in slice m except for the i-th spot, L mi Let be the number of neighbors of the i-th spot in slice m. Let Ψ be the conditional mean of the i-th spot neighbor in slice m. m ∈R q ×q Let be the conditional covariance matrix.

[0066] Step 3: Construct a dimensionality reduction model based on the gene expression matrix, initial dimensionality reduction parameters, and spatial dependencies between spots in the slice;

[0067] Step 4: Use an iterative optimization algorithm to estimate the parameters of the dimensionality reduction model and obtain the final low-dimensional representation of each slice.

[0068] Using the regularized gene expression matrix X m Taking modeling as an example, the dimensionality reduction model is specifically shown in equation (2):

[0069] x mi =α m +Wv mi +∈ mi (2)

[0070] Where, x mi ∈R p Let α be the gene expression vector of the i-th spot in the m-th slice. m ∈R pFor the sample-specific intercept term of the m-th slice, W∈R p×q The sample shared load matrix is ​​used to capture shared information between slices, v mi ∈R q A low-dimensional representation of the i-th spot in the m-th slice, used to capture biological information, ∈ mi ~N(0,Λ m ) represents the normally distributed random noise of the i-th spot in the m-th slice, where Λ m =diag(λ mi ,…,λ mp ) is a diagonal matrix.

[0071] Combining equations (1) and (2) yields the Gaussian version of the ProFAST model, with the parameters to be estimated including (W, α) m ,Λ m ,Ψ m (m∈{1,…,M}). Let X = (X1, ..., X) M ) T ∈R N×p V = (V1, ..., V) M ) T ∈R N×q ,in This is the total sample size. Assume there is a V... m The predicted value is defined as The pseudo-log-likelihood of the model can then be expressed as shown in equation (3):

[0072]

[0073] Wherein, P(x mi |v mi The result is obtained from equation (2). It is obtained from equation (1).

[0074] Gene expression matrix without regularization Taking modeling as an example, the dimensionality reduction model is specifically shown in equations (4) and (5):

[0075]

[0076]

[0077] in, Let a be the expression level of gene j in the i-th spot of slice m. mi The regularization factor for the i-th spot in slice m can be set to 1 or the sum of the counts of the i-th spot in slice m.j ∈R q Let α be the sample shared payload vector corresponding to gene j. mj f is the sample-specific intercept term corresponding to gene j in slice m. mij The Poisson rate of gene j in the i-th spot of slice m is an unknown quantity; v mi ∈R q Let m be the low-dimensional representation of the i-th spot in slice m, ∈ mij ~N(0,λ) mj ) represents the random noise of the normally distributed gene j in the i-th spot of slice m.

[0078] Let the load matrix Intercept vector α m =(α m1 ,…,α mp ), Combining equations (1), (4), and (5) yields the Poisson version of the ProFAST model, with the parameters to be estimated including (W, α) m ,Λ m ,Ψ m (m∈{1,…,M}). Suppose there exists a V m The predicted value is defined as The pseudo-log-likelihood of the model can then be expressed as shown in equation (6):

[0079]

[0080] in, From equation (4), lnP(f) is obtained. mij |v mi The result is obtained from equation (5). It is obtained from equation (1).

[0081] In practical applications, the ProFAST model can be iteratively optimized using the ICM-EM algorithm to update the model parameters. This algorithm alternates between ICM and EM steps until convergence, and finally outputs a low-dimensional representation of the slices, the model parameter estimates, and the likelihood function values.

[0082] The ProFAST model provided by this invention can simultaneously perform spatial dimensionality reduction on multiple slices. The model captures shared information between slices through a sample-shared payload matrix and captures biological information through sample-specific low-dimensional representations, enabling ProFAST to analyze multiple slices simultaneously. By establishing a continuous multivariate Hidden Markov Random Field on the low-dimensional representation, the spatial dependencies between spots in each slice are captured, allowing the model to achieve spatially aware dimensionality reduction.

[0083] The ProFAST model provided by this invention offers two versions of the algorithm: a Gaussian version and a Poisson version, which model the regularization matrix and the counting matrix, respectively. Furthermore, ProFAST boasts high computational efficiency. By using conditional autoregression to model local spatial dependencies, ProFAST's computational complexity is only O(n), making it suitable for multi-slice high-resolution SRT data.

[0084] Another embodiment of the present invention provides a dimensionality reduction system for spatial transcriptome data, the system comprising:

[0085] The data acquisition module is used to acquire the gene expression matrix and spatial coordinate matrix of multiple slices in the spatial transcriptome data;

[0086] The parameter determination module is used to determine the initial dimensionality reduction parameters and the spatial dependence between spots in each slice based on the gene expression matrix and the spatial coordinate matrix.

[0087] The model building module is used to construct a dimensionality reduction model based on the gene expression matrix, initial dimensionality reduction parameters, and spatial dependencies between spots in the slice.

[0088] The iterative optimization module is used to estimate the parameters of the dimensionality reduction model using an iterative optimization algorithm, and obtain the final low-dimensional representation of each slice.

[0089] For a detailed description of each module in the dimensionality reduction system described above, please refer to the description of each step in the dimensionality reduction method. It will not be repeated here. The dimensionality reduction system described above can achieve the same functionality as the dimensionality reduction method.

[0090] The above description is merely a few embodiments of this application and is not intended to limit this application in any way. Although this application discloses preferred embodiments as described above, it is not intended to limit this application. Any changes or modifications made by those skilled in the art without departing from the scope of the technical solution of this application using the disclosed technical content are equivalent to equivalent implementation cases and fall within the scope of the technical solution.

Claims

1. A method for dimensionality reduction of spatial transcriptome data, characterized in that, The method includes: Step 1, obtaining gene expression matrices and spatial coordinate matrices of multiple slices in spatial transcriptome data; Step 2, determining initial dimensionality reduction parameters and spatial dependencies between spots in each slice based on the gene expression matrices and spatial coordinate matrices; Step 3, constructing a dimensionality reduction model based on the gene expression matrices, the initial dimensionality reduction parameters, and the spatial dependencies between spots in each slice; Step 4, estimating the parameters of the dimensionality reduction model using an iterative optimization algorithm to obtain the final low-dimensional representation of each slice; Step 2 specifically includes: Step 21, determining initial dimensionality reduction parameters and the initial low-dimensional representation of each slice based on the gene expression matrices using principal component analysis; Step 22, determining the spatial dependencies between spots in each slice based on the spatial coordinate matrix and the initial low-dimensional representation of each slice; the spatial dependencies between spots in a slice are specifically defined as follows: given the low-dimensional representations of all spot samples other than the target spot in the slice, the low-dimensional representation of the target spot follows a normal distribution; wherein the mean of the normal distribution is the average of the low-dimensional representations of the neighbors of the target spot.

2. The method according to claim 1, characterized in that, Step 1 specifically includes: Step 11, obtaining the hypervariable genes and spatial coordinate matrix of each slice in the spatial transcriptome data to be analyzed, and integrating the hypervariable genes of each slice to obtain the integrated hypervariable genes; Step 12, obtaining the gene expression matrix of each slice based on the integrated hypervariable genes.

3. The method according to claim 2, characterized in that, Step 12 specifically includes: Step 121, extracting the counting matrix of the integrated hypervariable genes in each slice; Step 122, performing regularization on the counting matrix of each slice to obtain the gene expression matrix of each slice.

4. The method according to claim 3, characterized in that, The dimensionality reduction model is as follows: the gene expression matrix of each slice is the sum of three terms: the sample-specific intercept term of the slice, the product of the sample shared load matrix and the low-dimensional representation of the slice, and the random noise of the slice.

5. The method according to claim 2, characterized in that, The dimensionality reduction model is as follows: given the Poisson rate of the target spot in a slice, the gene expression of the target spot follows a Poisson distribution; and the logarithm of the Poisson rate of the target spot is the sum of three terms: the sample-specific intercept term of the slice, the product of the sample shared load vector and the low-dimensional representation of the slice, and the random noise of the slice.

6. The method according to claim 1, characterized in that, The iterative optimization algorithm is the ICM-EM algorithm.

7. The method according to claim 3, characterized in that, The regularization process is either logarithmic regularization or negative binomial regression regularization.

8. A dimensionality reduction system for spatial transcriptome data, characterized in that, The system includes: a data acquisition module for acquiring gene expression matrices and spatial coordinate matrices of multiple slices in spatial transcriptome data; a parameter determination module for determining initial dimensionality reduction parameters and spatial dependencies between spots in each slice based on the gene expression matrices and spatial coordinate matrices; a model construction module for constructing a dimensionality reduction model based on the gene expression matrices, the initial dimensionality reduction parameters, and the spatial dependencies between spots in each slice; and an iterative optimization module for estimating the parameters of the dimensionality reduction model using an iterative optimization algorithm to obtain the final low-dimensional representation of each slice. Specifically, the parameter determination module is used to: determine the initial dimensionality reduction parameters and the initial low-dimensional representation of each slice based on the gene expression matrices using principal component analysis; and determine the spatial dependencies between spots in each slice based on the spatial coordinate matrix and the initial low-dimensional representation of each slice. Specifically, the spatial dependencies between spots in a slice are defined as follows: given the low-dimensional representations of all spot samples other than the target spot in the slice, the low-dimensional representation of the target spot follows a normal distribution; wherein the mean of the normal distribution is the average of the low-dimensional representations of the target spot's neighbors.

Images