Method and device for deconvolution of spatial transcriptome data

Abstract

The invention discloses a spatial transcriptome data deconvolution method, which comprises the following steps of: (1) filtering spatial transcriptome data and single cell data, and reserving a part containing common gene expression information; (2) counting gene expression of each cell type based on the filtered single cell data and corresponding single cell annotation data, calculating the average gene expression level of each cell type, and constructing a reference matrix; (3) constructing an optimization model by taking space communication and space cell composition as decision variables; and (4) solving the optimization model by using an optimization algorithm with the aim of minimizing the difference between the predicted gene expression value and the real gene expression value to obtain joint inference of space deconvolution and a space communication mode. According to the method and the device, the deconvolution task of the space transcriptome data can be finished with high precision, and meanwhile, accurate inference of space communication is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of computational biology and bioinformatics, and in particular to a deconvolution method and device for spatial transcriptome data. Background Art

[0002] With the rapid development of high-throughput sequencing technology, people's research and understanding of gene expression have continued to deepen. The emergence of single-cell sequencing technology enables researchers to comprehensively capture the gene expression profile of individual cells and reveal their specific gene expression. However, due to the limitations of single-cell sequencing technology, this method is prone to losing the spatial location information of cells. This limitation has brought great trouble to scientific exploration when studying biological issues such as tumor microenvironment.

[0003] In order to make up for this deficiency, spatial transcriptomics technology came into being. Spatial transcriptomics can directly obtain the gene expression information of tissues and capture the gene expression profile in situ in tissues. This technology provides a new dimension for research while retaining spatial location information. However, limited by the resolution of spatial transcriptomics, the gene expression data obtained is usually a mixed result of multiple heterogeneous or homogeneous cell expressions. This mixed expression pattern poses challenges to analyzing the composition of tissue microenvironment and studying the development process and mechanism of cancer.

[0004] In order to analyze the cell composition in spatial transcriptome data, we need to complete the task of decomposing and inferring the mixed cells in each spatial location. This process requires clarifying the gene expression characteristics of each cell type, and single-cell sequencing technology can provide accurate expression data of individual cells. Therefore, combining single-cell data with spatial transcriptome data for spatial deconvolution analysis has become an important method for analyzing tissue microenvironment and exploring the mechanism of cancer development.

[0005] In previous studies, most researchers focused on the cell composition at each location in spatial transcriptome data, but ignored the rich location information contained in these data (Elosua-Bayes M, Nieto P, Mereu E, et al. SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes[J]. Nucleic Acids Research, 2021, 49 (9): e50-e50.). In the few deconvolution methods that consider location information, researchers usually fail to fully consider the potential interactions and communication relationships between cells (Dong R, Yuan GC. SpatialDWLS: accurate deconvolution of spatial transcriptomic data[J]. Genome Biology, 2021, 22 (1): 145.).

[0006] Therefore, to our knowledge, there is currently no method that can accurately model gene expression in spatial transcriptome data while achieving the dual tasks of deconvolution and spatial communication inference. This remains a technical challenge to be improved and overcome. Summary of the invention

[0007] The present invention provides a method and device for deconvolution of spatial transcriptome data, which can complete the deconvolution task of spatial transcriptome data with high precision and realize accurate inference of spatial communication at the same time.

[0008] The technical solution of the present invention is as follows: A deconvolution method for spatial transcriptome data, comprising: Step (1) filtering the spatial transcriptome data and single-cell data to retain the parts containing common gene expression information; Step (2) based on the filtered single-cell data and its corresponding single-cell annotation data, count the gene expression of each cell type, calculate the average gene expression level of each cell type, and construct a reference matrix; Step (3) uses spatial communication and spatial cell composition as decision variables to construct an optimization model; the constructed optimization model is as follows: Y = WY + XB + e In the formula, W and XThey are the space communication matrix and the cell composition matrix. B is the reference matrix, e is the difference between the predicted gene expression value and the actual gene expression; Y It is the real spatial transcriptome expression data; Step (4) aims to minimize the difference between the predicted gene expression value and the true gene expression value, and applies the optimization algorithm to solve the optimization model to obtain the joint inference of spatial deconvolution and spatial communication mode.

[0009] Step (1) includes: obtaining the intersection of genes in the spatial transcriptome data and genes in the single-cell data, and filtering the spatial transcriptome data and the single-cell data according to the intersection so that the two sets of data contain only common genes.

[0010] The gene filtering in step (1) is limited to retaining common genes, and highly variable genes are not screened or deleted. This strategy aims to preserve the information integrity of gene expression data to the greatest extent possible and avoid data loss or distortion due to excessive processing.

[0011] In order to complete the deconvolution task, it is necessary to construct reference data that can provide cell type expression characteristics. In step (2), the gene expression of each cell type is statistically analyzed, and the average gene expression level of each cell type is calculated. The obtained average expression value is used as the reference expression feature of the cell type. After obtaining the reference expression information of all cell types, this information is arranged in rows to form a matrix, which is the final reference matrix.

[0012] The spatial communication matrix and the cell composition matrix are both decision variables and need to be optimized and searched in the optimization model. In the optimization, the optimization direction is to minimize the difference e .

[0013] Step (4) includes: The difference is minimized with the constraints that the elements of the cell matrix are non-negative and the scope of spatial communication. e To optimize the direction, the optimization model is solved to obtain the joint inference of spatial deconvolution and spatial communication mode.

[0014] Preferably, in step (4), the coordinate information of the filtered spatial transcriptome data is used to calculate the spatial adjacency network D; and the spatial adjacency network D is used as a constraint condition in the optimization model.

[0015] Further, the calculation of the spatial adjacency network D using the coordinate information of the filtered spatial transcriptome data includes: Define each spatial point captured by sequencing as spot, for each element of the spatial adjacency network D dij Indicates j Is the spot the i Neighbor of a spot, d ij The value is 0 or 1, 0 represents non-neighbor relationship, 1 represents neighbor relationship; For each target spot, calculate the Euclidean distance between all other spots and the target spot, and set a cutoff threshold to screen whether they are adjacent.

[0016] Furthermore, the four other spots closest to the target spot are selected as neighbors.

[0017] In the present invention, both decision variables need to be optimized. However, in the optimization algorithm, optimizing multiple variables at the same time is prone to fall into the trouble of local optimal solution. Therefore, in the actual solution, the present invention adopts the sequential least squares algorithm (SLSQP) or the population differential evolution algorithm (DE) to solve the optimization model.

[0018] Preferably, solving the optimization model in step (4) includes: (i) Using regression model to optimize the initial cell composition matrix X 1 : Y = X 1 B + e (ii) Using initial cells to form a matrix X 1 Optimize space communication and obtain the initial space communication matrix W 1 : Y= W 1 Y+ X 1 B + e (iii) Using the initial value space communication matrix W 1 Continue to optimize the cell composition to obtain the final cell composition matrix X 2 : Y= W 1 Y+ X 2 B + e (iv) Completion of the cell matrix X 2 After solving the problem, bring it into the model to obtain the final spatial communication matrix W 2The solution: W 2 =(YX 2 B)Y -1 .

[0019] In the entire solution process of the optimization model, X 1. X 2. W 1 is obtained by optimization, and W 2 is obtained through the optimized inverse solution.

[0020] Preferably, the spatial transcriptome data deconvolution method of the present invention further comprises: verifying and evaluating the deconvolution result by using an evaluation index.

[0021] The evaluation index is the Pearson correlation coefficient between the spatially deconvolved cell composition vector a and the marker gene vector b. PCC , the calculation formula is as follows: .

[0022] The present invention also provides a deconvolution device for spatial transcriptome data, comprising: The data preprocessing module filters the spatial transcriptome data and single-cell data and retains the parts containing common gene expression information; The reference matrix construction module counts the gene expression of each cell type based on the filtered single-cell data and its corresponding single-cell annotation data, calculates the average gene expression level of each cell type, and constructs a reference matrix; The optimization model construction module uses space communication and space cell composition as decision variables to construct an optimization model; the constructed optimization model is as follows: Y = WY + XB + e In the formula, W and X They are the space communication matrix and the cell composition matrix. B is the reference matrix, e is the difference between the predicted gene expression value and the actual gene expression; Y It is the real spatial transcriptome expression data; The optimization model solving module aims to minimize the difference between the predicted gene expression value and the true gene expression value, and applies the optimization algorithm to solve the optimization model to obtain the joint inference of spatial deconvolution and spatial communication mode.

[0023] Compared with the prior art, the present invention has the following beneficial effects: The present invention proposes a new deconvolution method for spatial transcriptome data, which combines spatial communication inference and cell composition inference, so that the method of the present invention can complete two tasks at the same time. In addition, the deconvolution method of spatial transcriptome data of the present invention fully considers spatial information, overcoming the problem that the existing deconvolution method ignores spatial information. In addition, the present invention also models spatial communication information, quantifies spatial communication, and improves the accuracy of deconvolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 Schematic diagram of the process of the deconvolution calculation method of spatial transcriptome data of the present invention. DETAILED DESCRIPTION

[0025] The present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be pointed out that the embodiments described below are intended to facilitate the understanding of the present invention and do not have any limiting effect on the present invention.

[0026] It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical terms and scientific terms) used herein have the same meaning as generally understood by those skilled in the art in the field to which the present invention belongs. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with the meanings in the context of the prior art, and will not be interpreted with idealized or overly formal meanings unless defined as herein.

[0027] The numbers of the steps mentioned in each embodiment are only for the convenience of description and have no substantial sequence relationship. Different steps in each specific embodiment can be combined in different sequences to achieve the purpose of the invention.

[0028] like Figure 1 As shown, a spatial deconvolution method based on least squares optimization specifically includes the following steps: 101) Preprocessing steps of spatial transcriptome data and single-cell data: To preprocess all input data, the present invention requires four types of data: spatial transcriptome data, spatial transcriptome coordinate information, single-cell data, and single-cell annotation data. Spatial transcriptome data is gene expression data obtained through spatial transcriptome technology. Spatial transcriptome coordinate information is the coordinate data of the spatial transcriptome capture position. Single-cell data is single-cell data that matches the spatial transcriptome organization. Single-cell annotation data is the cell type annotation data for all cells in the single-cell sequencing data. In the preprocessing step, we obtained the intersection of the genes in the spatial transcriptome data and the genes in the single-cell data, and filtered the two sets of data with this intersection so that only common genes are contained in the data. Here, this method does not select highly variable genes to retain the original characteristics of the data to the greatest extent.

[0029] 102) Reference matrix construction step: Use single-cell data and single-cell annotation data to calculate the reference matrix. Specifically, in this step, this method will count the expression of each cell type in the single-cell data and calculate the average expression of each cell type. The calculated average expression is used as an approximation of the expression of this cell type. Then, the expression of these cell types is arranged in rows to construct the reference matrix. B After the calculation is completed, B The rows of the matrix represent cell types and the columns represent genes.

[0030] 103) Construction and solution steps of the optimization model: The standard form of the optimization model consists of three parts, namely constraints, decision variables and optimization objectives. In this step, there are two parts of constraints, which come from two characteristics of biological analysis. Since the analysis goal of the deconvolution method is to analyze the composition of cell types in space, and the cell composition in the tissue cannot be negative, the first constraint is the non-negative constraint of the cell composition matrix. Secondly, since the biological tissues of multicellular organisms have spatial continuity, cell communication and information transmission are more likely to occur between spatial positions that are close in space. Taking this feature into account, the second part of the constraint of this model is the constraint on the scope of communication within the tissue, and only considers the communication of local cells. Specifically, this method will calculate the spatial adjacency network D , and restrict such a neighbor network to spatial communication. In spatial transcriptome data, we call each spatial point captured by sequencing a spot. Spatial neighbor network D Describes the adjacency relationship of all spots, so the spatial adjacency network D is a square matrix, the number of rows and columns of which are the number of spots. D Every element of d ij , d ij Indicates whether the jth spot is a neighbor of the ith spot, and its value is 0 or 1. A value of 0 represents a non-neighbor relationship, and a value of 1 represents a neighbor relationship. When calculating the spatial adjacency network, for each spot, we calculate the Euclidean distance between all other spots and this spot, and set a truncation threshold to screen whether they are adjacent. In this method, we select the four other spots closest to the target spot as neighbors, that is, in the spatial adjacency network D In each row vector of , only four elements are 1 and the other elements are 0. Complete the spatial adjacency network DAfter calculating, this method will use the Hadamard product to set constraints on spatial communication. In the model assumption, spots will affect each other's expression. This algorithm models this influence information into a linear model. In actual expression, the effect of cells on their own expression should be the maximum value, so the absolute value of the linear influence coefficient should also be limited to less than 1 in the linear model.

[0031] Different from the past methods, which only focus on deconvolution tasks, thus ignoring the rich spatial information and information transfer, the model is as follows: Y = XB + e In the formula, B is the reference matrix constructed in step 102), X Represents the space cell composition matrix, e Represents the difference between the gene expression predicted by the model and the actual gene expression.

[0032] Furthermore, the model was improved to model spatial information and cell composition. The specific optimization model is as follows: Y = WY + XB + e In the formula, Y is the real spatial transcriptome expression data, which is an expression matrix. The rows represent spots, the columns represent genes, and the matrix element (i, j) represents the expression level of gene j in the i-th spot. B is the reference matrix constructed in step 102), W and X is the decision variable of the model, where W Represents the spatial communication strength, X Represents the space cell composition matrix, e Represents the difference between the gene expression predicted by the model and the actual gene expression. This model will be optimized multiple times, first using the initial solution of the initial model solution, then modeling spatial communication, optimizing communication and cell composition. In the optimization model, we modeled communication expression and gene expression, and used a linear model to calculate the predicted value of expression, that is, WY + XB Therefore, based on such modeling, the optimization goal of this optimization problem is to optimize the predicted value and the true expression value, that is, to minimize the residual e .

[0033] Furthermore, this model will use the Sequential Least Squares Programming (SLSQP) algorithm or the Differential Evolution Algorithm (DE) algorithm to solve the optimization model according to the settings. The SLSQP algorithm converts the original optimization problem into its quadratic subproblem, performs a first-order approximation on the constraints, and converts the original optimization problem into a quadratic programming problem for convex optimization, which can effectively avoid the influence of local optimality. In addition, the DE algorithm is a heuristic algorithm that searches for the optimal solution by simulating population evolution. This optimization method can effectively avoid the influence of local optimal values. The details are as follows: Y = X 1 B + e First, the regression model is used to optimize the initial cell composition X 1. Then, using the initial composition X 1. Optimize the initial value of space communication W 1: Y= W 1 Y+ X 1 B + e After completing the initial inference of space communication, bring in the initial inference W 1. Continue to optimize the cell composition to obtain the final cell composition X 2: Y= W 1 Y+ X 2 B + e After solving the cell composition, it is brought into the model for the final spatial communication inference: W 2 =(YX 2 B)Y -1 During the entire optimization process of the model, X 1, X 2, W 1 is obtained by optimization, and W 2 is obtained through the optimized inverse solution.

[0034] Through these optimization methods, this method will obtain the cell composition matrix X Communication Matrix with Space W .

[0035] 104) Result verification step: After completing the inference of the cell composition of the spatial transcriptome data, this method provides a verification method for optimizing the results. Specifically, the evaluation index is the Pearson correlation coefficient between the spatial transcriptome cell composition results and the marker genes, calculated as follows:

[0036] In the formula, a , b are two vectors, the distribution vector of cell types and the expression vector relative to marker genes. you and sb The vectors are a With vector b The distribution vector of cell types is derived from the cell composition results. X The marker gene vector b It was obtained from spatial transcriptome data.

[0037] Although the above methods are illustrated and described as a series of actions for simplicity of explanation, it should be understood and appreciated that these methods are not limited by the order of the actions, because according to one or more embodiments, some actions may occur in a different order and / or concurrently with other actions from those illustrated and described herein or not illustrated and described herein but understandable to those skilled in the art.

[0038] Those skilled in the art will further appreciate that the various illustrative logic blocks, modules, circuits, and algorithm steps described in conjunction with the embodiments disclosed herein may be implemented as electronic hardware, computer software, or a combination of the two. To clearly illustrate this interchangeability of hardware and software, various illustrative components, blocks, modules, circuits, and steps are generally described above in terms of their functionality. Whether such functionality is implemented as hardware or software depends on the specific application and the design constraints imposed on the overall system. The technician may implement the described functionality in different ways for each specific application, but such implementation decisions should not be interpreted as resulting in a departure from the scope of the present invention.

[0039] The various illustrative logic blocks, modules, and circuits described in conjunction with the embodiments disclosed herein may be implemented or executed with a general purpose processor, a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA) or other programmable logic device, discrete gate or transistor logic, discrete hardware components, or any combination thereof designed to perform the functions described herein. A general purpose processor may be a microprocessor, but in an alternative, the processor may be any conventional processor, controller, microcontroller, or state machine. The processor may also be implemented as a combination of computing devices, such as a combination of a DSP and a microprocessor, a plurality of microprocessors, one or more microprocessors in cooperation with a DSP core, or any other such configuration.

[0040] The steps of the method or algorithm described in conjunction with the embodiments disclosed herein may be embodied directly in hardware, in a software module executed by a processor, or in a combination of the two. The software module may reside in a RAM memory, a flash memory, a ROM memory, an EPROM memory, an EEPROM memory, a register, a hard disk, a removable disk, a CD-ROM, or any other form of storage medium known in the art. An exemplary storage medium is coupled to a processor so that the processor can read and write information from / to the storage medium. In an alternative, a storage medium may be integrated into a processor. The processor and the storage medium may reside in an ASIC. The ASIC may reside in a user terminal. In an alternative, the processor and the storage medium may reside in a user terminal as discrete components.

[0041] In one or more exemplary embodiments, the functions described may be implemented in hardware, software, firmware, or any combination thereof. If implemented as a computer program product in software, each function may be stored on or transmitted by a computer-readable medium as one or more instructions or codes. Computer-readable media include both computer storage media and communication media, including any media that facilitates the transfer of a computer program from one place to another. Storage media may be any available media that can be accessed by a computer. As an example and not limitation, such a computer-readable medium may include RAM, ROM, EEPROM, CD-ROM or other optical disk storage, disk storage or other magnetic storage device, or any other medium that can be used to carry or store the desired program code in the form of an instruction or data structure and can be accessed by a computer. Any connection is also properly referred to as a computer-readable medium. For example, if the software is transmitted from a website, server, or other remote source using a coaxial cable, fiber optic cable, twisted pair, digital subscriber line (DSL), or wireless technologies such as infrared, radio, and microwaves, the coaxial cable, fiber optic cable, twisted pair, DSL, or wireless technologies such as infrared, radio, and microwaves are included in the definition of the medium. Disk and disc as used herein include compact disc (CD), laser disc, optical disc, digital versatile disc (DVD), floppy disk and Blu-ray disc, where disks usually reproduce data magnetically, while discs reproduce data optically with lasers. Combinations of the above should also be included within the scope of computer-readable media.

[0042] The embodiments described above provide a detailed description of the technical solutions and beneficial effects of the present invention. It should be understood that the above are only specific embodiments of the present invention and are not intended to limit the present invention. Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A deconvolution method for spatial transcriptome data, characterized in that: include: Step (1) filtering the spatial transcriptome data and single-cell data to retain the parts containing common gene expression information; Step (2) based on the filtered single-cell data and its corresponding single-cell annotation data, count the gene expression of each cell type, calculate the average gene expression level of each cell type, and construct a reference matrix; Step (3) uses spatial communication and spatial cell composition as decision variables to construct an optimization model; the constructed optimization model is as follows: Y = WY + XB + ε In the formula, W and X They are the space communication matrix and the cell composition matrix. B is the reference matrix, ε is the difference between the predicted gene expression value and the actual gene expression; Y It is the real spatial transcriptome data; Step (4) aims to minimize the difference between the predicted gene expression value and the true gene expression value, and applies the optimization algorithm to solve the optimization model to obtain the joint inference of spatial deconvolution and spatial communication mode.

2. The deconvolution method of spatial transcriptome data according to claim 1, characterized in that Step (1) includes: obtaining the intersection of genes in the spatial transcriptome data and genes in the single-cell data, and filtering the spatial transcriptome data and the single-cell data according to the intersection so that the two sets of data contain only common genes.

3. The deconvolution method of spatial transcriptome data according to claim 1, characterized in that Step (4) includes: The difference is minimized with the constraints that the elements of the cell matrix are non-negative and the scope of spatial communication. ε To optimize the direction, the optimization model is solved to obtain the joint inference of spatial deconvolution and spatial communication mode.

4. The deconvolution method of spatial transcriptome data according to claim 3, characterized in that: In step (4), the coordinate information of the filtered spatial transcriptome data is used to calculate the spatial adjacency network D; the spatial adjacency network D is used as a constraint condition in the optimization model.

5. The deconvolution method of spatial transcriptome data according to claim 4, characterized in that: The calculation of the spatial adjacency network D using the coordinate information of the filtered spatial transcriptome data includes: Define each spatial point captured by sequencing as spot, for each element of the spatial adjacency network D d ij Indicates j Is the spot the i Neighbor of a spot, d ij The value is 0 or 1, 0 represents non-neighbor relationship, 1 represents neighbor relationship; For each target spot, calculate the Euclidean distance between all other spots and the target spot, and set a cutoff threshold to screen whether they are adjacent.

6. The deconvolution method of spatial transcriptome data according to claim 5, characterized in that: Select the four other spots closest to the target spot as neighbors.

7. The deconvolution method of spatial transcriptome data according to claim 1, characterized in that: Solving the optimization model in step (4) includes: (i) Using regression model to optimize the initial cell composition matrix X 1: Y = X 1 B + ε (ii) Using initial cells to form a matrix X 1 Optimize space communication and obtain the initial space communication matrix W 1: Y = W 1 Y+X 1 B + ε (iii) Using the initial value space communication matrix W 1 Continue to optimize the cell composition to obtain the final cell composition matrix X 2: Y = W 1 Y+X 2 B + ε (iv) Completion of the cell matrix X After solving 2, it is brought into the model to obtain the final space communication matrix W Solution for 2: W 2 =(YX 2 B)Y -1 。 8. The method for deconvolution of spatial transcriptome data according to claim 1, characterized in that: Also includes: The deconvolution results are verified and evaluated through evaluation indicators.

9. A deconvolution device for spatial transcriptome data, characterized in that: include: The data preprocessing module filters the spatial transcriptome data and single-cell data and retains the parts containing common gene expression information; The reference matrix construction module counts the gene expression of each cell type based on the filtered single-cell data and its corresponding single-cell annotation data, calculates the average gene expression level of each cell type, and constructs a reference matrix; The optimization model construction module uses space communication and space cell composition as decision variables to construct an optimization model; the constructed optimization model is as follows: Y = WY + XB + ε In the formula, W and X They are the space communication matrix and the cell composition matrix. B is the reference matrix, ε is the difference between the predicted gene expression value and the actual gene expression; Y It is the real spatial transcriptome data; The optimization model solving module takes minimizing the difference between the predicted gene expression value and the true gene expression value as the objective function, applies the optimization algorithm to solve the optimization model, and obtains the joint inference of spatial deconvolution and spatial communication mode.

Images