Cell region determination method and device, electronic equipment and storage medium

Abstract

The embodiment of the invention provides a cell region determination method and device, electronic equipment and a storage medium, and belongs to the technical field of cell segmentation. The cell region determination method comprises the following steps: acquiring gene expression quantity information from spatial transcription data of a target section; wherein the gene expression quantity information comprises position data of a plurality of space sequencing sites and the gene expression quantity of each space sequencing site; calculating a corresponding target distribution parameter when the gene expression quantity of each spatial sequencing site obeys preset distribution; calculating a first probability that each spatial sequencing site belongs to a cell coverage area according to the target distribution parameters of adjacent sites of each spatial sequencing site; and determining a cell region of the target section according to the first probability and the position data. According to the embodiment of the invention, the accuracy of determining the cell region can be improved.

Description

Technical Field

[0001] The present disclosure relates to the technical field of cell segmentation, and in particular to a method and device for determining a cell region, an electronic device, and a storage medium. Background Art

[0002] At present, determining the cell region is an important step in conducting spatial cell research. In the related art, methods for determining cell regions include methods for determining cell regions based on staining images obtained based on H&E and other technologies, and methods for delineating cells through partial RNA signals. Among them, the method based on staining images is based only on image features, without considering RNA signals, resulting in a low average RNA capture amount in cells obtained after cell segmentation, which is not conducive to the downstream analysis of spatial omics data. The method for delineating cells based on partial RNA signals requires a large amount of a certain RNA signal, that is, the method depends on the capture amount of a specific gene, so the method is not suitable for sequencing-based spatial omics data with sparse RNA signals. Therefore, there is currently a lack of a method that can accurately determine the cell region and thus achieve reasonable cell segmentation. Summary of the invention

[0003] The embodiments of the present disclosure provide a method and device for determining a cell region, an electronic device, and a storage medium, which can improve the accuracy of determining a cell region.

[0004] To achieve the above object, a first aspect of an embodiment of the present application provides a method for determining a cell region, the method comprising:

[0005] Obtaining gene expression information from the spatial transcription data of the target slice; wherein the gene expression information includes position data of multiple spatial sequencing sites and the gene expression level of each of the spatial sequencing sites;

[0006] Calculating the target distribution parameters corresponding to when the gene expression amount of each of the spatial sequencing sites follows a preset distribution;

[0007] Calculating a first probability that each of the spatial sequencing sites belongs to a cell coverage area according to the target distribution parameters of the adjacent sites of each spatial sequencing site; wherein the adjacent sites of the spatial sequencing sites are spatial sequencing sites adjacent to the spatial sequencing sites;

[0008] The cell region of the target slice is determined according to the first probability and the position data.

[0009] In some embodiments, a first log-likelihood function of the preset distribution is constructed according to the gene expression level of each of the spatial sequencing sites;

[0010] Determine the hidden variable according to the first log-likelihood function, and execute the target step cyclically, wherein the target step includes: calculating the conditional probability of the hidden variable according to the randomly initialized original distribution parameters; maximizing the first log-likelihood function according to the conditional probability to obtain preliminary distribution parameters; if the target step is executed a preset number of times, use the preliminary distribution parameters when the preset number of times is reached as the target distribution parameters, or if the preliminary distribution parameters converge, use the preliminary distribution parameters as the target distribution parameters.

[0011] In some embodiments, constructing the first log-likelihood function of the preset distribution according to the gene expression level of each of the spatial sequencing sites includes:

[0012] Constructing the probability expression data of the preset distribution according to the gene expression amount of each of the spatial sequencing sites;

[0013] Constructing a state variable of the spatial sequencing site according to the probabilistic expression data; wherein the state variable is used to indicate whether the spatial sequencing site belongs to a cell coverage area or a non-cell coverage area;

[0014] Constructing a distribution variable of the preset distribution according to the probability expression data; wherein the distribution variable is used to indicate that the type of the preset distribution is a cell coverage area type or a non-cell coverage area type;

[0015] The first log-likelihood function is constructed according to the state variable and the distribution variable.

[0016] In some embodiments, the conditional probability includes a cell region probability, a first variable probability, and a second variable probability, and the original distribution parameters include a first distribution parameter and a second distribution parameter;

[0017] The step of calculating the conditional probability of the latent variable according to the randomly initialized original distribution parameters includes:

[0018] Calculate the probability that the spatial sequencing site belongs to the cell coverage area according to the preset conditions and the original distribution parameters to obtain the cell area probability;

[0019] Perform logarithmic calculation according to the first distribution parameter to obtain the first variable probability;

[0020] Perform logarithmic differentiation calculation according to the second distribution parameter to obtain the second variable probability.

[0021] In some embodiments, maximizing the first log-likelihood function according to the conditional probability to obtain preliminary distribution parameters includes:

[0022] Summing the cell region probabilities to obtain a probability sum;

[0023] Calculating the product sum of the cell region probability and the second variable probability, performing ratio calculation on the product sum and the accumulated value of the cell region probability to obtain a first probability ratio;

[0024] A second probability ratio is calculated according to the cell region probability, the first variable probability, and the second variable probability;

[0025] An expectation function of the first log-likelihood function is calculated, and the expectation function is maximized according to the probability sum, the first probability ratio, and the second probability ratio to obtain the preliminary distribution parameters.

[0026] In some embodiments, determining the target distribution parameter corresponding to when the gene expression level of each of the spatial sequencing sites follows a preset distribution includes:

[0027] Constructing a second probability density function of the preset distribution according to the gene expression level of each of the spatial sequencing sites;

[0028] Obtaining observation data of the spatial sequencing site belonging to the cell coverage area in a preset number of times, and constructing a second likelihood function according to the observation data and the second probability density function;

[0029] Performing logarithmic processing on the second likelihood function to obtain a logarithmic likelihood function;

[0030] The log-likelihood function is derivatized to obtain the target distribution parameters.

[0031] In some embodiments, determining the first probability that the spatial sequencing site belongs to the cell coverage area according to the target distribution parameter of the adjacent site includes:

[0032] Calculate the adjacent message data between the adjacent site and the spatial sequencing site according to the target distribution parameter;

[0033] A first probability that the spatial sequencing site belongs to a cell coverage area is determined based on the neighboring message data.

[0034] In some embodiments, the calculating the adjacent message data between the adjacent site and the spatial sequencing site according to the target distribution parameter includes:

[0035] Calculating the conditional probability that the spatial sequencing site belongs to the cell coverage area according to the target distribution parameters of the spatial sequencing site;

[0036] Calculating a first potential energy difference of the adjacent message according to the first state data of the adjacent site, the target distribution parameter of the adjacent site, and a preset potential function; wherein the first state data indicates whether the adjacent site belongs to a cell coverage area or a non-cell coverage area;

[0037] The adjacent message data is calculated according to the conditional probability and the first potential energy difference.

[0038] In some embodiments, determining the first probability that the spatial sequencing site belongs to a cell coverage area according to the neighboring message data includes:

[0039] Summing a plurality of adjacent message data to obtain total message data;

[0040] Calculating a second potential energy difference of the spatial sequencing site according to the second state data of the spatial sequencing site and the target distribution parameter of the spatial sequencing site; wherein the second state data indicates that the spatial sequencing site belongs to a cell coverage area or a non-cell coverage area;

[0041] A first probability that the spatial sequencing site belongs to a cell coverage area is calculated based on the total message data and the second potential energy difference.

[0042] In some embodiments, determining the cell region of the target slice according to the first probability and the position data includes:

[0043] constructing a cell mask matrix according to the first probability and the position data, and constructing a spatial count matrix according to the gene expression information;

[0044] The cell area of ​​the target slice is determined according to the cell mask matrix and the spatial count matrix.

[0045] In some embodiments, constructing a cell mask matrix according to the first probability and the position data comprises:

[0046] Determining a second probability that the spatial sequencing site belongs to a non-cell coverage area according to the target distribution parameters of the adjacent sites;

[0047] The first probability is compared with the second probability, and the cell mask matrix is ​​constructed according to the comparison result, the preset mask value, and the position data.

[0048] In some embodiments, the root determines the cell region of the target slice according to the cell mask matrix and the spatial count matrix, comprising:

[0049] Constructing a convolution kernel according to a preset cell diameter, and performing convolution calculation according to the convolution kernel and the spatial counting matrix to obtain a convolution matrix;

[0050] Performing dot product calculation on the convolution matrix and the cell mask matrix to obtain a cell label matrix; wherein the row and column dimensions of the cell label matrix are the same as those of the spatial count matrix, the cell label matrix includes the position data of the spatial sequencing sites and the cell labels corresponding to the spatial sequencing sites, and the cell labels include cell region labels;

[0051] The cell region of the target slice is determined according to the position data corresponding to the cell region label.

[0052] To achieve the above-mentioned purpose, a second aspect of an embodiment of the present application provides a cell region determination device, the device comprising:

[0053] An information acquisition unit, used to acquire gene expression information from the spatial transcription data of the target slice; wherein the gene expression information includes position data of multiple spatial sequencing sites and the gene expression level of each of the spatial sequencing sites;

[0054] A distribution parameter calculation unit, used to calculate the target distribution parameter corresponding to the gene expression amount of each of the spatial sequencing sites when it follows a preset distribution;

[0055] A probability calculation unit, used to calculate a first probability that each of the spatial sequencing sites belongs to a cell coverage area according to the target distribution parameters of the adjacent sites of each spatial sequencing site; wherein the adjacent sites of the spatial sequencing sites are spatial sequencing sites adjacent to the spatial sequencing sites;

[0056] A cell region determining unit is used to determine the cell region of the target slice according to the first probability and the position data.

[0057] To achieve the above-mentioned purpose, a third aspect of an embodiment of the present application proposes an electronic device, including a memory and a processor, wherein the memory stores a computer program, and when the processor executes the computer program, the method described in the first aspect is implemented.

[0058] To achieve the above-mentioned purpose, the fourth aspect of the embodiments of the present application proposes a computer-readable storage medium, which stores a computer program. When the computer program is executed by a processor, the method described in the first aspect is implemented.

[0059] To achieve the above-mentioned purpose, the fifth aspect of the embodiments of the present application proposes a computer program product, which includes a computer program. The computer program is read and executed by a processor of a computer device, so that the computer device executes the method described in the first aspect.

[0060] The cell region determination method and device, electronic device and storage medium proposed in the embodiment of the present application, the target distribution parameter of the preset distribution is calculated according to the gene expression amount of each spatial sequencing site, and the first probability that the spatial sequencing site belongs to the cell coverage area is calculated according to the target distribution parameter of the adjacent site of the spatial sequencing site, so as to determine the cell region of the target slice according to the first probability and the position data of the spatial sequencing site. It can be seen from this that the example of the present application is to determine the cell region of the target slice completely based on the gene expression amount. Compared with the method based on the staining image in the related art, the gene expression information obtained by the embodiment of the present application directly to the spatial transcription data is processed, and the situation that the RNA signal amount of the demarcated cell region is low is reduced, and the consistency of the RNA signal and the cell region determination result is enhanced. Compared with the method for determining the cell region according to partial RNA, the embodiment of the present application can reduce the dependence on specific genes and improve the universality of the cell region determination method. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] The accompanying drawings are used to provide further understanding of the technical solution of the present disclosure and constitute a part of the specification. Together with the embodiments of the present disclosure, they are used to explain the technical solution of the present disclosure and do not constitute a limitation on the technical solution of the present disclosure.

[0062] Figure 1 is a flow chart of a method for determining a cell region provided in an embodiment of the present application;

[0063] Figure 2 yes Figure 1 Flow chart of step S102 in FIG.

[0064] Figure 3 yes Figure 2 Flow chart of step S201 in FIG.

[0065] Figure 4 yes Figure 2 Flow chart of calculating hidden variable conditional probability in step S202;

[0066] Figure 5 yes Figure 2 Flow chart of calculating preliminary distribution parameters in step S202;

[0067] Figure 6 yes Figure 1 A flowchart of another embodiment of step S102 in FIG.

[0068] Figure 7 yes Figure 1 Flow chart of step S103 in FIG.

[0069] Figure 8 yes Figure 7 Flow chart of step S701 in FIG.

[0070] Fig. 9 yes Figure 7 Flow chart of step S702 in FIG.

[0071] Fig.10 yes Figure 1 Flow chart of step S104 in FIG.

[0072] Fig.11 yes Fig.10 Flow chart of step S1001 in FIG.

[0073] Fig.12 yes Fig.10 Flowchart of step S1002 in FIG.

[0074] Fig.13 It is a schematic diagram of the cell area determined by the method of artificial labeling based on RNA signals;

[0075] Fig.14 is a schematic diagram of a cell region determined using the method provided in the examples of the present application;

[0076] Fig.15 It is a schematic diagram of the cell area determined using an image-based reference method (such as the Stardist method);

[0077] Fig.16 It is a schematic diagram of the cell region determined by the method of delineating the cell region based on partial RNA signals (such as the Baysor method);

[0078] Fig.17 is a schematic diagram of a cell region determination device provided in an embodiment of the present application;

[0079] Fig.18 It is a schematic diagram of the hardware structure of the electronic device provided in the embodiment of the present application. DETAILED DESCRIPTION

[0080] In order to make the purpose, technical solution and advantages of the present disclosure more clear, the present disclosure is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present disclosure and are not used to limit the present disclosure.

[0081] Before further describing the embodiments of the present disclosure in detail, the nouns and terms involved in the embodiments of the present disclosure are described. The nouns and terms involved in the embodiments of the present disclosure are subject to the following interpretations:

[0082] Spatial Transcriptomics (ST) sequencing technology: a technology that analyzes RNA from a spatial level, used to analyze all RNA in a single tissue section. There are four capture areas on each slide used for library construction by 10x Genomics spatial transcriptomics, and the size of each capture area is 6.5x6.5mm. Each capture area contains 5,000 barcoded spots, each with a diameter of 55μm and a center-to-center distance of 100μm. Each spot includes multiple capture probes that can bind to RNA, and each probe carries unique spatial barcodes to mark the spatial position of the captured RNA. Through on-machine sequencing, the sequence of each RNA transcript can be mapped back to the original position in the tissue section, thereby providing data support for downstream tasks. It can be seen that ST sequencing technology is a technology that can simultaneously obtain RNA expression and RNA spatial location information in a single experiment.

[0083] Spatial Transcriptomics with Error-Corrected Barcoding and Sequencing (Stereo-seq): including the design of DNA nanoball (DNB) patterned array chips, in situ sequencing to determine the spatial coordinates of unique barcode oligonucleotides, point connection of UMI-polyT oligonucleotides, in situ capture of RNA from tissues, cDNA amplification, library construction sequencing and data analysis. Among them, DNB sequencing is based on a patterned array of in situ sequencing, which can serve as the basis for spatial resolution transcription technology with high resolution and large field of view. Compared with HDST (HDST is a high-resolution spatial transcriptome technology that carves a large number of 2.05μm holes on the surface of the slide and distributes silica magnetic beads in the holes. The magnetic beads are used to capture the mRNA of the cells in the corresponding position for reverse transcription, library construction and transcriptome sequencing), Slide-seqV2 (Slide-seqV2 is a high-sensitivity near-cellular resolution spatial transcriptome analysis that can achieve better RNA capture efficiency through steps such as magnetic bead synthesis and array indexing), Visium (Visium is a spatial transcriptome sequencing technology that integrates gene expression with immunohistochemical images of tissue sections to locate the gene expression levels of different cells in the tissue to the original spatial position of the tissue), DBiT-seq (DBiT-seq is a spatial multi-omics sequencing technology based on microfluidic barcode labeling), etc., Stereo-seq has a high resolution of 100μm. 2With more spots, smaller spot size and center-to-center distance, Stereo-seq has higher resolution.

[0084] At present, for spatial transcriptome sequencing technology with subcellular resolution, determining the cell region of tissue sections is an important step in conducting spatial cell research. In related technologies, methods for determining cell regions include the following two categories: the first category is to determine the cell region based on staining images obtained by technologies such as H&E; the second category is to delineate the cell region based on partial RNA.

[0085] Among them, the methods based on stained images include the following two: the first is to identify the stained images based on computer vision feature extraction, random forest and other machine learning models (such as the Ilastik algorithm, which is an image segmentation and analysis method based on machine learning) to determine the cell area. The second is to identify the stained images based on deep learning models of convolutional neural networks (such as the DeepCell model, the Stardist model (Stardist model is a model for object detection and segmentation), etc.) to determine the cell area.

[0086] Methods for delineating cell regions based on partial RNA (such as the Baysor algorithm (the Baysor algorithm is an algorithm based on the MRF segmentation idea) and the ClusterMap algorithm (the ClusterMap algorithm is a clustering-based algorithm)) use image-based in situ transcription sequencing methods such as MERFISH to stain a small number of specific genes, and then delineate cells based on the staining images, that is, determine the cell region.

[0087] Among the above methods, the method based on machine learning model can automatically detect cell regions in fluorescent staining images, but it requires a labeled data set for model training. In this case, the trained machine learning model may have poor generalization ability to other data sets, that is, the accuracy of cell region recognition is low when applied to other data sets.

[0088] The method based on deep learning model requires manual adjustment of the results to obtain better cell area determination results, and the determined cell area is different from the actual cell distribution. This is because the image needs to be collected through subsequent experiments processed by ST, which may cause the staining image and RNA signal to be inconsistent.

[0089] The method of delineating cell regions based on partial RNA requires a high capture amount of specific genes if accurate cell delineation results are expected. Therefore, this method is not suitable for processing sequencing-based spatial transcriptome data with sparse RNA signals. In addition, this method takes a long time to run and consumes a lot of memory.

[0090] Based on this, the embodiments of the present application propose a method and device for determining a cell region, an electronic device, and a storage medium, which can improve the accuracy of determining a cell region.

[0091] The cell region determination method provided in the embodiments of the present application is described below.

[0092] Reference Figure 1 In some embodiments, the cell region determination method provided in the embodiments of the present application includes but is not limited to steps S101 to S104.

[0093] Step S101, obtaining gene expression information from the spatial transcription data of the target slice; wherein the gene expression information includes position data of multiple spatial sequencing sites and the gene expression level of each spatial sequencing site;

[0094] Step S102, calculating the target distribution parameters corresponding to when the gene expression level of each spatial sequencing site follows a preset distribution;

[0095] Step S103, calculating a first probability that each spatial sequencing site belongs to a cell coverage area according to a target distribution parameter of an adjacent site of each spatial sequencing site; wherein an adjacent site of a spatial sequencing site is a spatial sequencing site adjacent to the spatial sequencing site;

[0096] Step S104: determining a cell region of the target slice according to the first probability and the position data.

[0097] Steps S101 to S104 shown in the embodiment of the present application are calculated according to the target distribution parameters of the preset distribution of each spatial sequencing site, and the first probability that the spatial sequencing site belongs to the cell coverage area is calculated according to the target distribution parameters of the adjacent sites of the spatial sequencing site, so as to determine the cell area of ​​the target slice according to the first probability and the position data of the spatial sequencing site. It can be seen from this that the present application example is to determine the cell area of ​​the target slice completely based on the gene expression amount. Compared with the method based on the staining image in the related art, the gene expression information obtained by the spatial transcription data is directly processed by the embodiment of the present application, and the situation that the RNA signal amount of the demarcated cell area is low is reduced, and the consistency of the RNA signal and the cell area determination result is enhanced. Compared with the method for determining the cell area according to partial RNA, the embodiment of the present application can reduce the dependence on specific genes and improve the universality of the cell area determination method.

[0098] In step S101 of some embodiments, the target slice refers to a frozen embedded tissue slice to be subjected to ST treatment to determine the cell region. The target slice is a biological slice obtained in accordance with relevant legal provisions, and the biological slice includes a plant slice and an animal slice. The target slice is subjected to ST treatment to obtain spatial transcription data, and gene expression information is obtained according to the spatial transcriptome data. Gene expression information is used to describe the position data of the spatial sequencing site, and the gene expression amount corresponding to the spatial sequencing site. Among them, the spatial sequencing site refers to the spot in the ST technology, the position information is used to indicate the position of the spot in the corresponding capture area, and the gene expression amount refers to the aggregate value of the RNA expression amount captured by the spot from the corresponding position of the target slice. For example, after ST treatment, a list can be obtained, and the list is used to describe the sub-expression amounts of various genes captured by each spatial sequencing site. The sub-expression amounts of various genes are aggregated (i.e., summed) to obtain the gene expression amount of the corresponding spatial sequencing site. It can be understood that the gene expression information can be represented in the form of a matrix, the corresponding position information can be represented in the form of (x, y) coordinates, and the gene expression amount is the value of the (x, y) coordinate in the matrix.

[0099] In step S102 of some embodiments, the preset distribution refers to a pre-set probability distribution, and the preset distribution can be a negative binomial distribution, a Poisson distribution, a Poisson distribution with a zero expansion term, etc. For ease of explanation, in the following embodiments, the preset distribution is taken as an example of a negative binomial distribution. The negative binomial distribution is used to describe the relationship between success and failure in an independent event. For example, in the embodiments of the present application, whether each spatial sequencing site belongs to a cell coverage area can be regarded as an independent event. According to the negative binomial distribution, it can be determined that each spatial sequencing site belongs to a cell coverage area (i.e., success) and belongs to a non-cell coverage area (i.e., failure), and then the probability of each spatial sequencing site belonging to a cell coverage area can be estimated. It should be understood that when the preset distribution is set to other probability distributions and the adaptive modifications made to the embodiments of the present application should also belong to the protection scope of the present application. The target distribution parameter refers to the distribution parameter corresponding to the preset distribution when the gene expression amount of each spatial sequencing site obeys the preset distribution. It can be understood that the gene expression amount of each spatial sequencing site can be subject to an independent preset distribution, so the target distribution parameter corresponding to the preset distribution can be calculated according to the gene expression amount.

[0100] Reference Figure 2 In some embodiments, step S102 includes but is not limited to steps S201 to S202.

[0101] Step S201, constructing a first log-likelihood function of a preset distribution according to the gene expression level of each spatial sequencing site;

[0102] Step S202, determining the hidden variable according to the first log-likelihood function, and executing the target step in a loop, the target step including: calculating the conditional probability of the hidden variable according to the randomly initialized original distribution parameters; maximizing the first log-likelihood function according to the conditional probability to obtain the preliminary distribution parameters; if the target step is executed a preset number of times, the preliminary distribution parameters when the preset number of times is reached are used as the target distribution parameters, or if the preliminary distribution parameters converge, the preliminary distribution parameters are used as the target distribution parameters.

[0103] In step S201 of some embodiments, the first log-likelihood function is a function of a parameter of a preset distribution. The first log-likelihood function for each spatial sequencing site can be constructed based on the gene expression level of the spatial sequencing site.

[0104] Reference Figure 3 In some embodiments, step S201 includes but is not limited to steps S301 to S304.

[0105] Step S301, constructing probability expression data of a preset distribution according to the gene expression amount of each spatial sequencing site;

[0106] Step S302, constructing a state variable of the spatial sequencing site according to the probabilistic expression data; wherein the state variable is used to indicate whether the spatial sequencing site belongs to a cell coverage area or a non-cell coverage area;

[0107] Step S303, constructing a distribution variable of a preset distribution according to the probability expression data; wherein the distribution variable is used to indicate that the type of the preset distribution is a cell coverage area type or a non-cell coverage area type;

[0108] Step S304: construct a first log-likelihood function according to the state variables and the distribution variables.

[0109] In step S301 of some embodiments, when the preset distribution is a negative binomial distribution as an example, for each spatial sequencing site, an expression (i.e., probabilistic expression data) as shown in the following formula (1) can be constructed based on the gene expression amount.

[0110] P(X|p 0 ,p 1 ,r,r′,θ,θ′)=p 1 NB(X|r,θ)+p 0 ·NB(X|r′,θ′)...Equation (1)

[0111] Among them, X represents the gene expression level of the spatial sequencing site, p 1 Indicates the probability that the gene expression level belongs to the foreground signal expression level, p 0Indicates the probability that the gene expression level belongs to the background signal gene expression level. The foreground signal refers to the RNA signal corresponding to the cell area, and the background signal refers to the RNA signal due to the technical noise parameter. Technical noise can include transcript positioning error, gene expression estimation error, etc. Transcript positioning error can be caused by factors such as uncertainty in the measurement method and changes in signal intensity. Gene expression estimation error can be caused by factors such as noise in the RNA signal, limitations of data processing and analysis methods, etc. r and θ represent the parameters of the preset distribution corresponding to the foreground signal. r' and θ' represent the parameters of the preset distribution corresponding to the background signal.

[0112] In steps S302 to S304 of some embodiments, in order to facilitate calculation and improve the calculation speed of the EM algorithm, dummy variables are constructed based on the expression shown in formula (1) to separate the variables in the expression. Specifically, state variables S, distribution variables M, variables Y and variables Z are constructed, wherein variables Y and variables Z are proposed for the convenience of calculation and have no practical significance. The state variable S is used to represent the state of the spatial sequencing site, that is, it indicates that the spatial sequencing site belongs to the cell coverage area or the non-cell coverage area. The distribution variable M represents the type of preset distribution, that is, it indicates which set of parameters the spatial sequencing site obeys. The cell coverage area type represents the preset distribution of the set of parameters (r, θ), and the non-cell area type represents the preset distribution of the set of parameters (r′, θ′). The first log-likelihood function shown in the following formula (2) is constructed based on the state variable S, distribution variable M, variable Y and variable Z.

[0113]

[0114] Where T represents the set of state variables S. When the spatial sequencing site includes multiple states, it can be represented by the set T. k represents the kth group of distribution parameters, including two groups of distribution parameters (r,θ) and (r′,θ′).

[0115] The value of variable Y can be calculated according to the following formula (3), and the value of variable Z can be calculated according to the following formula (4).

[0116]

[0117]

[0118] The benefit of step S301 to step S303 is that, since the log-likelihood function of the classic negative binomial distribution contains a quadratic gamma function, directly using the EM algorithm to calculate the parameters of such a negative binomial distribution will result in the problem of being unable to provide an accurate numerical solution, and the EM algorithm calculation will consume a lot of computing time. The embodiment of the present application solves the above problem by introducing dummy variables, thereby improving the calculation speed of the subsequent EM algorithm.

[0119] In step S202 of some embodiments, the target distribution parameters can be calculated according to the EM algorithm. Among them, the EM algorithm includes an E step (i.e., an expectation step) and an M step (i.e., a maximization step). By looping the target step, the E step and the M step are performed alternately until convergence, and the target distribution parameters are obtained. It can be understood that the target step can be regarded as a process of performing an E step and an M step. In the target step, the hidden variable refers to a variable that appears in the first log-likelihood function but cannot be directly observed. In the E step, the conditional probability of the hidden variable is calculated according to the parameter value of the randomly initialized preset distribution (i.e., the original distribution parameter). In the M step, the first log-likelihood function is maximized based on the conditional probability calculated in the E step to update the parameters of the preset distribution, that is, to obtain a new parameter value of the preset distribution (i.e., the initial distribution parameter). It can be understood that the distribution parameters used in the next round of the E step are the calculation results of the previous round of the M step (i.e., the initial distribution parameters). The goal of the M step is to find the parameter value that maximizes the first log-likelihood function. Repeat the above target steps until convergence. Convergence means that the number of repeated executions reaches a preset number of times, or the initial distribution parameters converge (i.e., the initial distribution parameters are stable and no longer change significantly). The corresponding preliminary distribution parameters at the time of convergence are used as target distribution parameters. It is understandable that the specific values ​​of the preset number of times and the randomly initialized original distribution parameters can be adaptively set according to actual needs, and this embodiment of the present application does not make specific limitations on this.

[0120] Reference Figure 4 In some embodiments, the conditional probability includes the cell region probability, the first variable probability and the second variable probability, and the original distribution parameters include the first distribution parameters and the second distribution parameters. In step S202, "calculating the conditional probability of the hidden variable according to the randomly initialized original distribution parameters" includes but is not limited to steps S401 to S403.

[0121] Step S401, calculating the probability that the spatial sequencing site belongs to the cell coverage area according to the preset conditions and the original distribution parameters, and obtaining the cell area probability;

[0122] Step S402, performing logarithmic calculation according to the first distribution parameter to obtain a first variable probability;

[0123] Step S403, performing logarithmic differentiation calculation according to the second distribution parameter to obtain the second variable probability.

[0124] It should be noted that step S401 to step S403 can be regarded as step E in the EM algorithm. The specific calculation process of step E is described below.

[0125] In step S401 of some embodiments, the cell region probability refers to the probability that the spatial sequencing site belongs to the cell coverage area. The cell region probability is calculated according to the preset conditions, the original distribution parameters and the following formula (5):

[0126] The preset condition refers to a certain value (i.e. x t ) and a set of references (i.e. ) is a restriction. h represents the iteration of each step. For example, the original distribution parameter randomly initialized when the E step is performed for the first time is represented by θ (0) 、r (0) After the original distribution parameters are updated according to the initial distribution parameters obtained in the M step, the distribution parameters for the second E step are expressed as θ (1) 、r (1) . represents the first distribution parameter, represents the second distribution parameter.

[0127] In some embodiments, in step S402 to step S403, the first variable probability and the second variable probability have no practical meaning and are for the convenience of calculating the iterative data introduced for the variables Y and Z. The first variable probability is calculated according to the first distribution parameter and the following formula (6): The second variable probability is calculated based on the second distribution parameter and the following formula (7):

[0128]

[0129]

[0130] in:

[0131]

[0132]

[0133] Formula (8) represents the differential of the logarithm of the γ function, and formula (9) represents the conversion relationship between λ, θ, and r.

[0134] Reference Figure 5In some embodiments, in step S202, "maximizing the first log-likelihood function according to the conditional probability to obtain preliminary distribution parameters" includes but is not limited to steps S501 to S504.

[0135] Step S501, summing up the cell region probabilities to obtain a probability sum;

[0136] Step S502, calculating the product sum of the cell region probability and the second variable probability, and performing ratio calculation on the product sum and the accumulated value of the cell region probability to obtain a first probability ratio;

[0137] Step S503, calculating a second probability ratio according to the cell region probability, the first variable probability, and the second variable probability;

[0138] Step S504, calculating the expected function of the first log-likelihood function, maximizing the expected function according to the probability sum, the first probability ratio, and the second probability ratio to obtain preliminary distribution parameters.

[0139] It should be noted that step S501 to step S504 can be regarded as the M step in the EM algorithm. The specific calculation process of the M step is described below.

[0140] In steps S501 to S503 of some embodiments, the cell region probability is calculated according to the following formula (10): Perform the summation to obtain the probability and The cell region probability is calculated according to the following formula (11): With the second variable probability The product sum of the cell area probability The first probability ratio is obtained by calculating the accumulated value of According to the probability of cell area First variable probability Second variable probability The second probability ratio is calculated as follows:

[0141]

[0142]

[0143]

[0144] In step S504 of some embodiments, the expectation of the first log-likelihood function is calculated with reference to the following formula (13) to obtain an expected function (i.e., Q function). The probability sum, the first probability ratio, and the second probability ratio are substituted into the expected function, and the expected function is maximized according to the gradient direction of the expected function to obtain the parameters corresponding to the maximization of the expected function (i.e., preliminary distribution parameters).

[0145] Q(φ,φ h )=E[log(S,M,Y,Z;φ)|x,φ (h) ]=Q(p,p (h) )+Q(φ,φ (h) )......Formula (13)

[0146] Among them, Q(p,p (h) ) represents the Q function for the distribution variable M, Q(φ,φ (h) ) represents the Q function for the variable of the parameter of the preset distribution. Q(p,p (h) ) The specific expression is shown in the following formula (14), Q(φ,φ (h) )The specific expression is shown in formula (15).

[0147]

[0148]

[0149]

[0150] Reference Figure 6 In some other embodiments, step S102 also includes but is not limited to steps S601 to S604.

[0151] Step S601, constructing a second probability density function of a preset distribution according to the gene expression level of each spatial sequencing site;

[0152] Step S602, obtaining observation data of the spatial sequencing sites belonging to the cell coverage area in a preset number of times, and constructing a second likelihood function according to the observation data and the second probability density function;

[0153] Step S603, performing logarithmic processing on the second likelihood function to obtain a logarithmic likelihood function;

[0154] Step S604, performing a derivative process on the log-likelihood function to obtain target distribution parameters.

[0155] It should be noted that in step S601 to step S604, the parameters of the preset distribution (i.e., the target distribution parameters) are solved according to the maximum likelihood probability. It is understandable that in addition to the EM algorithm and the maximum likelihood probability, other methods such as the maximum a posteriori probability can also be used to solve the target distribution parameters. The specific calculation process of the maximum likelihood probability is described below.

[0156] In step S601 of some embodiments, a corresponding second probability density function is constructed according to the expression of the negative binomial distribution and the gene expression level of each spatial sequencing site.

[0157] In step S602 of some embodiments, the observation data refers to the result that the spatial sequencing site belongs to the cell coverage area in the preset number of experiments. A result set can be constructed based on the observation data. For example, assuming that the preset number of times is 10 times, the result set can be expressed as {1, 0, 1, 1, 0, 1, 0, 1, 1}, where 1 represents the result that the spatial sequencing site belongs to the cell coverage area, and 0 represents the result that the spatial sequencing site belongs to the non-cell coverage area. A second likelihood function is constructed based on the result set and the second probability density function.

[0158] In step S603 of some embodiments, in order to simplify calculation, the logarithm of the second likelihood function is taken to obtain a logarithmic likelihood function.

[0159] In step S604 of some embodiments, the log-likelihood function is maximized to solve the target distribution parameter of the preset distribution. The method for maximizing the log-likelihood function may be to find the derivative of the log-likelihood function. The target distribution parameter is the parameter value that makes the log-likelihood function reach the maximum value.

[0160] In step S103 of some embodiments, the adjacent site refers to the spot adjacent to the spatial sequencing site, and the adjacent site refers to the distance between the adjacent site and the spatial sequencing site within a preset distance range, and the distance can be calculated based on the position data. The specific value of the preset distance range can be adaptively set according to actual needs, for example, it can be adaptively set according to the accuracy of determining the cell area, and the embodiments of the present application are not specifically limited to this. The first probability that the corresponding spatial sequencing site belongs to the cell coverage area is calculated based on the target distribution parameters of the adjacent site. It is understandable that belonging to the cell coverage area refers to the cell area of ​​the target slice corresponding to the spatial sequencing site.

[0161] Reference Figure 7 In some embodiments, step S103 includes but is not limited to steps S701 to S702.

[0162] Step S701, calculating the adjacent message data between the adjacent sites and the spatial sequencing sites according to the target distribution parameters;

[0163] Step S702: determining a first probability that the spatial sequencing site belongs to a cell coverage area according to the adjacent message data.

[0164] It should be noted that when the target distribution parameters (including r, r′, θ, θ′) are calculated according to step S102, the conditional probability that each spatial sequencing site belongs to the cell coverage area can actually be calculated according to the following formula (16), thereby determining some high-confidence sites (i.e., spatial sequencing sites with large conditional probabilities). However, in actual situations, it is unlikely that most areas of a cell correspond to high-confidence sites. If the cell area is determined directly based on these high-confidence sites, the obtained cell area is incomplete. Therefore, it is necessary to use a method (such as the method shown in steps S701 to S702) to use the sites in the high-confidence site concentration area as sites of the cell coverage area so as to obtain a complete cell area.

[0165]

[0166] The conditional probability calculated according to formula (16) can be used to construct a confidence matrix. Each spatial sequencing site can be represented as a point in the confidence matrix. The position data of the spatial sequencing site can be represented as the coordinate data of the corresponding point. The gene expression level of the spatial sequencing site can be represented as the corresponding value of the corresponding point in the confidence matrix.

[0167] It can be understood that the method shown in step S701 to step S702 is actually a belief propagation algorithm (BP algorithm). The basis of the BP algorithm is based on the Markov Random Field (MRF), which is used to describe the correlation of all nodes in the field (i.e., spatial sequencing sites). In the BP algorithm, the joint probability distribution of each spatial sequencing site and the adjacent site can be assumed to be as shown in the following formula (17).

[0168]

[0169] Among them, x represents the spatial sequencing site, y represents the adjacent site of the spatial sequencing site, and x (i,j) Indicates the gene expression level at the spatial sequencing site, y (i,j) Indicates the status of the spatial sequencing site, y (k,l) represents the state of the adjacent sites, and N represents the set of adjacent sites.

[0170] Each spatial sequencing site has a state value (such as whether it belongs to the state of the cell coverage area) and an observed value (such as the gene expression of the spatial sequencing site), and the likelihood function between the state value and the observed value reflects the dependency between the two. The potential energy between adjacent spatial sequencing sites can reflect the correlation between adjacent spatial sequencing sites. The BP algorithm uses messages between adjacent spatial sequencing sites to update the labeling state of the entire MRF. The BP algorithm is an iterative algorithm. After multiple iterations, if the confidence of all spatial sequencing sites (i.e., the first probability of belonging to the cell coverage area) no longer changes, the state of each spatial sequencing site at this time is regarded as the optimal state. Steps S701 to S702 are described below.

[0171] In step S701 of some embodiments, a spatial sequencing site is randomly selected, the adjacent site of the spatial sequencing site is determined according to the above method, and the adjacent message data is calculated according to the target distribution parameter of the adjacent site. The adjacent message data refers to the message transmitted by the adjacent site to the spatial sequencing site, and the message includes information such as what distribution and state the spatial sequencing site should be in.

[0172] Reference Figure 8 In some embodiments, step S701 includes but is not limited to steps S801 to S803.

[0173] Step S801, calculating the conditional probability that the spatial sequencing site belongs to the cell coverage area according to the target distribution parameters of the spatial sequencing site;

[0174] Step S802, calculating a first potential energy difference of the adjacent message according to the first state data of the adjacent site, the target distribution parameter of the adjacent site, and a preset potential function; wherein the first state data indicates whether the adjacent site belongs to a cell coverage area or a non-cell coverage area;

[0175] Step S803, calculating adjacent message data according to the conditional probability and the first potential energy difference.

[0176] In step S801 of some embodiments, the method for calculating the conditional probability that the spatial sequencing site belongs to the cell coverage area according to the target distribution parameter is referred to the above formula (16), X (i,j) Represents the gene expression level at the spatial sequencing site.

[0177] In step S802 of some embodiments, the first potential energy difference represents the difference between the adjacent messages belonging to high confidence and low confidence. The first potential energy difference can be calculated according to the following formula (18).

[0178] λ (k,l) =log(ψ′(x (k,l) ,1) / ψ′(x (k,l),0))......Formula (18)

[0179] Among them, ψ′(x (k,l) ,1) indicates that the adjacent message is of high confidence, ψ′(x (k,l) ,0) indicates that the adjacent message is of low confidence. (k,l) ,1) and ψ′(x (k,l) ,0) can be calculated according to the potential energy equation (i.e., potential function) shown in equation (19).

[0180]

[0181] When ψ′(x (k,l) ,1) and ψ′(x (k,l) ,0), replace y in formula (19) (i,j) is replaced by the first state data of the adjacent site, and r, θ, r′, θ′ are replaced by the target distribution parameters of the adjacent site. The first state data is used to indicate whether the adjacent site belongs to the cell-covered area or the non-cell-covered area. It can be understood that in formula (19), y (i,j) =1 means it belongs to the cell coverage area, y (i,j) =0 indicates the area is not covered by cells.

[0182] In step S803 of some embodiments, the adjacent message data l is calculated according to the conditional probability, the first potential energy difference and the following formula (20): (k,l)→(i,j) .

[0183]

[0184] Among them, p represents the conditional probability, and (o,p) represents other adjacent sites of the spatial sequencing site (i,j) except the adjacent site (k,l).

[0185] It can be understood that the adjacent message data shown in formula (20) is actually the message data m shown in the following formula: (k,l)→(i,j) The log-likelihood form of is:

[0186]

[0187] In addition, the marginal probability distribution can be expressed as shown below:

[0188]

[0189] In step S702 of some embodiments, a first probability that the spatial sequencing site belongs to a cell coverage area is determined based on the neighbor message data.

[0190] It is understandable that steps S701 to S702 represent only one iteration, and it is necessary to select another spatial sequencing site and repeat steps S701 to S702 until the first probabilities of all spatial sequencing sites no longer change. The first probability when no longer changes is used as the confidence that each spatial sequencing site ultimately belongs to the cell coverage area.

[0191] Reference Fig. 9 In some embodiments, step S702 includes but is not limited to steps S901 to S903.

[0192] Step S901, summing up multiple adjacent message data to obtain total message data;

[0193] Step S902, calculating a second potential energy difference of the spatial sequencing site according to the second state data of the spatial sequencing site and the target distribution parameter of the spatial sequencing site; wherein the second state data indicates whether the spatial sequencing site belongs to a cell coverage area or a non-cell coverage area;

[0194] Step S903, calculating a first probability that the spatial sequencing site belongs to the cell coverage area according to the total message data and the second potential energy difference.

[0195] In steps S901 to S903 of some embodiments, the adjacent message data of multiple adjacent sites corresponding to the spatial sequencing site are summed to obtain the total message data. Referring to the method shown in step S801 to step S803, the second potential energy difference λ of the spatial sequencing site is calculated. (i,j) The first probability of the spatial sequencing site is calculated according to the total message data, the second potential energy difference and the following formula (21). The first probability represents the probability that the spatial sequencing site belongs to the cell coverage area.

[0196]

[0197] The advantage of steps S801 to S803 and steps S901 to S903 is that the message data m (k,l)→(i,j) It can be seen from the calculation formula that the message data m (k,l)→(i,j)It is actually the result of the aggregation of the adjacent message data of all adjacent sites of a spatial sequencing site. Wherein, the adjacent message data refers to the data that the adjacent site should be in what kind of distribution, what kind of state, etc., which is passed to the spatial sequencing site. It can be seen that the embodiment of the present application is based on the state of the adjacent site to estimate the state of each spatial sequencing site. Therefore, for the spatial sequencing site that is directly determined as a non-high confidence position according to the conditional probability, the final state of the spatial sequencing site can be further re-estimated according to the adjacent site of the spatial sequencing site. In this way, when there are multiple high-confidence spatial sequencing sites in the adjacent site of the spatial sequencing site, the spatial sequencing site can also be determined as a high-confidence spatial sequencing site. That is, the spatial sequencing site of the edge position of the high-confidence region can also be included in the high-confidence region to avoid omission. At this time, a complete cell region can be obtained according to multiple spatial sequencing sites belonging to the cell coverage area.

[0198] In step S104 of some embodiments, multiple spatial sequencing sites can be divided according to the first probability to obtain high-confidence sites belonging to the cell coverage area. The enclosed area of ​​the high-confidence site can be determined based on the position data of the high-confidence site, and the enclosed area is the area in the capture area corresponding to the target slice cell area, so the cell area of ​​the target slice can be determined.

[0199] Reference Fig.10 In some embodiments, step S104 includes but is not limited to steps S1001 to S1002.

[0200] Step S1001, constructing a cell mask matrix according to the first probability and position data, and constructing a spatial count matrix according to gene expression information;

[0201] Step S1002, determining the cell region of the target slice according to the cell mask matrix and the spatial count matrix.

[0202] In step S1001 of some embodiments, a cell mask matrix M is constructed according to the first probability and position data of each spatial sequencing site. Specifically, it is determined whether the spatial sequencing site belongs to the cell coverage area according to the first probability. If it belongs to the cell coverage area, the position data of the spatial sequencing site in the cell mask matrix M can be assigned a value of 1. If it belongs to the non-cell coverage area, the position data of the spatial sequencing site in the cell mask matrix M can be assigned a value of 0. A spatial count matrix X is constructed according to the position data of each spatial sequencing site and the gene expression amount of the spatial sequencing site. In the spatial count matrix X, the gene expression amount is used as the assignment of the corresponding position data.

[0203] Reference Fig.11 In some embodiments, step S1001 includes but is not limited to steps S1101 to S1102.

[0204] Step S1101, determining a second probability that the spatial sequencing site belongs to a non-cell coverage area according to target distribution parameters of adjacent sites;

[0205] Step S1102, comparing the first probability with the second probability, and constructing a cell mask matrix according to the comparison result, the preset mask value, and the position data.

[0206] In step S1101 of some embodiments, the second probability that the spatial sequencing site belongs to the non-cell coverage area is calculated based on the target distribution parameters of the adjacent sites and the following formula (22). It can be understood that the specific method for calculating the second probability can refer to the method for calculating the first probability, and this embodiment of the application will not be repeated.

[0207]

[0208] In step S1102 of some embodiments, the first probability P(y (i,j) =1) and the second probability P(y (i,j) =0) for comparison, and determine the assignment of the position data of each spatial sequencing site in the cell mask matrix M according to the comparison result. Specifically, if the first probability is greater than the second probability, it indicates that the corresponding spatial sequencing site belongs to the cell coverage area, and the position data of the spatial sequencing site is assigned to 1. If the first probability is less than the second probability, it indicates that the corresponding spatial sequencing site belongs to the non-cell coverage area, and the position data of the spatial sequencing site is assigned to 0.

[0209] In step S1002 of some embodiments, cell segmentation is performed based on the cell mask matrix and the spatial count matrix to ultimately determine multiple cell regions of the target slice.

[0210] Reference Fig.12 In some embodiments, step S1002 includes but is not limited to steps S1201 to S1203.

[0211] Step S1201, constructing a convolution kernel according to a preset cell diameter, and performing convolution calculation according to the convolution kernel and the spatial counting matrix to obtain a convolution matrix;

[0212] Step S1202, performing dot product calculation on the convolution matrix and the cell mask matrix to obtain a cell label matrix; wherein the row and column dimensions of the cell label matrix are the same as those of the spatial count matrix, the cell label matrix includes the position data of the spatial sequencing sites and the cell labels corresponding to the spatial sequencing sites, and the cell labels include cell region labels;

[0213] Step S1203, determining the cell region of the target slice according to the position data corresponding to the cell region label.

[0214] In step S1201 of some embodiments, a convolution kernel is constructed according to a preset cell diameter. The cell diameter can be determined according to prior knowledge, such as determining the cell diameter corresponding to each cell in the target slice according to prior knowledge, and then calculating the average of multiple cell diameters. Convolution kernel G is constructed according to the average of multiple cell diameters. Convolution kernel G can be an odd convolution kernel. For example, if the cell diameter s is 15, then convolution kernel G=s*s. It can be understood that when the image is represented, the unit of cell diameter is pixel point. The convolution kernel G and the spatial count matrix X are convolved as shown in the following formula (23) to obtain the convolution matrix V.

[0215] V=-Conv(X,G)......Formula (23)

[0216] In step S1202 of some embodiments, the convolution matrix V is dot-producted with the cell mask matrix M to obtain a cell label matrix L. The row and column dimensions of the cell label matrix L are the same as those of the spot in the ST process, that is, when the gene expression information is in matrix form, the row and column dimensions of the cell label matrix L are the same as the gene expression information. The cell label matrix L is used to describe whether the position data corresponding to each spatial sequencing site belongs to the cell coverage area. When it belongs to the cell coverage area, the cell label corresponding to the spatial sequencing site is the cell area label. When it belongs to the non-cell coverage area, the cell label corresponding to the spatial sequencing site is the non-cell area label. It is understandable that the cell area label can be in the form of a number, and a number of 0 represents a non-cell area label, and a number greater than 0 represents a cell area label. In the cell area label, different number values ​​represent corresponding different cell areas.

[0217] In step S1203 of some embodiments, the cell region of the target slice is obtained by segmentation according to the position data of the cell region labels in the cell label matrix L.

[0218] The advantage of steps S1201 to S1203 is that not only can the cell region of the target slice be determined, but also each cell region can be distinguished according to the number corresponding to the cell region label.

[0219] In a specific embodiment, referring to Figures 13 to 16 Based on the cell region determination method provided in the embodiments of the present application, the cell region is determined for the mouse brain spatial transcriptome data produced by the Stereo-seq sequencing technology. Figures 13 to 16 This is a diagram showing the effect of determining the cell area in a certain area captured in a real experiment. Fig.13 It is a cell region determined by artificial labeling of RNA signals. Fig.14 is a cell region determined using the method provided in the embodiment of the present application, Fig.15is the cell area determined using image-based reference methods (such as the Stardist method), Fig.16 It is a cell region determined by a method that delineates cell regions based on partial RNA signals (such as the Baysor method). Figures 13 to 16 By comparing, Fig.16 The method shown has the disadvantage of segmenting one cell region into multiple cell regions, i.e. Fig.16 The method shown in the present application embodiment is not a complete cell region. And the method shown in the present application embodiment is close to the result of the artificial labeling method, and the scope of the method cell delineation shown in the present application embodiment is large, there is no discrete scattered cell, and the adjacent area can be integrated into a cell region. In addition, the cell delineated by the method shown in the present application embodiment is closer to the size of the real cell in theory. It can be seen from this that compared to the method for determining the cell region based on image or based on partial RNA, the method provided in the present application embodiment is more suitable for analyzing high-throughput spatial transcriptome sequencing data, and the cell region determination result reliability obtained by the method provided in the present application embodiment is high, and the RNA signal of the delineated cell region is strong, which is conducive to downstream data analysis.

[0220] The cell region determination device provided in the embodiment of the present application is described below.

[0221] Reference Fig.17 In some embodiments, the present application also provides a cell region determination device, the cell region determination device comprising:

[0222] The information acquisition unit 1701 is used to acquire gene expression information from the spatial transcription data of the target slice; wherein the gene expression information includes the position data of multiple spatial sequencing sites and the gene expression amount of each spatial sequencing site;

[0223] A distribution parameter calculation unit 1702 is used to calculate the target distribution parameter corresponding to the gene expression amount of each spatial sequencing site when it follows a preset distribution;

[0224] A probability calculation unit 1703 is used to calculate a first probability that each spatial sequencing site belongs to a cell coverage area according to a target distribution parameter of an adjacent site of each spatial sequencing site;

[0225] The cell region determining unit 1704 is used to determine the cell region of the target slice according to the first probability and the position data.

[0226] It can be seen that the contents of the above-mentioned cell region determination method embodiment are all applicable to the embodiment of the present cell region determination device. The functions specifically implemented by the present cell region determination device embodiment are the same as those in the above-mentioned cell region determination method embodiment, and the beneficial effects achieved are also the same as those achieved by the above-mentioned cell region determination method embodiment.

[0227] Reference Fig.18 , Fig.18 The hardware structure of an electronic device of another embodiment is illustrated, and the electronic device includes:

[0228] The processor 1801 may be implemented by a general-purpose CPU (Central Processing Unit), a microprocessor, an application-specific integrated circuit (Application Specific Integrated Circuit, ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of the present application;

[0229] The memory 1802 can be implemented in the form of a read-only memory (ROM), a static storage device, a dynamic storage device, or a random access memory (RAM). The memory 1802 can store an operating system and other application programs. When the technical solution provided in the embodiments of this specification is implemented by software or firmware, the relevant program code is stored in the memory 1802, and the processor 1801 calls and executes the cell region determination method of the embodiment of this application;

[0230] Input / output interface 1803, used to implement information input and output;

[0231] The communication interface 1804 is used to realize the communication interaction between the device and other devices. The communication can be realized through a wired manner (such as USB, network cable, etc.) or a wireless manner (such as mobile network, WIFI, Bluetooth, etc.);

[0232] A bus 1805 that transmits information between the various components of the device (e.g., the processor 1801, the memory 1802, the input / output interface 1803, and the communication interface 1804);

[0233] The processor 1801 , the memory 1802 , the input / output interface 1803 and the communication interface 1804 are connected to each other in communication within the device via the bus 1805 .

[0234] The embodiment of the present application further provides a computer program product, which includes a computer program. A processor of a computer device reads and executes the computer program, so that the computer device executes and implements the above-mentioned cell region determination method.

[0235] The terms "first", "second", "third", "fourth", etc. (if any) in the specification of the present disclosure and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence. It should be understood that the data used in this way can be interchangeable where appropriate, so that the embodiments of the present disclosure described herein can, for example, be implemented in an order other than those illustrated or described herein. In addition, the terms "comprises" and "comprising" and any variations thereof are intended to cover non-exclusive inclusions, for example, a process, method, system, product or device that includes a series of steps or units is not necessarily limited to those steps or units that are clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0236] It should be understood that in the present disclosure, "at least one (item)" means one or more, and "plurality" means two or more. "And / or" is used to describe the association relationship of associated objects, indicating that three relationships may exist. For example, "A and / or B" can mean: only A exists, only B exists, and A and B exist at the same time, where A and B can be singular or plural. The character " / " generally indicates that the objects associated before and after are in an "or" relationship. "At least one of the following" or similar expressions refers to any combination of these items, including any combination of single or plural items. For example, at least one of a, b or c can mean: a, b, c, "a and b", "a and c", "b and c", or "a and b and c", where a, b, c can be single or multiple.

[0237] It should be understood that in the description of the embodiments of the present application, the meaning of multiple (or multiple items) is more than two, greater than, less than, exceed, etc. are understood to not include the number, and above, below, within, etc. are understood to include the number.

[0238] In the several embodiments provided in the present disclosure, it should be understood that the disclosed systems, devices and methods can be implemented in other ways. For example, the device embodiments described above are only schematic. For example, the division of units is only a logical function division. There may be other division methods in actual implementation. For example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be an indirect coupling or communication connection through some interfaces, devices or units, which can be electrical, mechanical or other forms.

[0239] The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed on multiple network units. Some or all of the units may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0240] In addition, each functional unit in each embodiment of the present disclosure may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit. The above-mentioned integrated unit may be implemented in the form of hardware or in the form of software functional units.

[0241] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present disclosure is essentially or the part that contributes to the prior art or all or part of the technical solution can be embodied in the form of a software product, and the computer software product is stored in a storage medium, including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the various embodiments of the present disclosure. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (Read-Only Memory, referred to as ROM), random access memory (Random Access Memory, referred to as RAM), disk or optical disk and other media that can store program codes.

[0242] It should also be understood that the various implementations provided in the embodiments of the present application can be combined arbitrarily to achieve different technical effects.

[0243] The above is a specific description of the implementation methods of the present disclosure, but the present disclosure is not limited to the above implementation methods. Technical personnel familiar with the art can also make various equivalent modifications or substitutions without violating the spirit of the present disclosure. These equivalent modifications or substitutions are all included in the scope defined by the claims of the present disclosure.

Claims

1. A method for determining a cell region, characterized in that: The method comprises: Obtaining gene expression information from the spatial transcription data of the target slice; wherein the gene expression information includes position data of multiple spatial sequencing sites and the gene expression level of each of the spatial sequencing sites; Calculating the target distribution parameters corresponding to when the gene expression amount of each of the spatial sequencing sites follows a preset distribution; Calculating a first probability that each of the spatial sequencing sites belongs to a cell coverage area according to the target distribution parameters of the adjacent sites of each spatial sequencing site; wherein the adjacent sites of the spatial sequencing sites are spatial sequencing sites adjacent to the spatial sequencing sites; The cell region of the target slice is determined according to the first probability and the position data.

2. The method according to claim 1, characterized in that The calculating of the target distribution parameters corresponding to when the gene expression amount of each of the spatial sequencing sites obeys a preset distribution includes: Constructing a first log-likelihood function of the preset distribution according to the gene expression level of each of the spatial sequencing sites; Determine the hidden variable according to the first log-likelihood function, and execute the target step cyclically, wherein the target step includes: calculating the conditional probability of the hidden variable according to the randomly initialized original distribution parameters; maximizing the first log-likelihood function according to the conditional probability to obtain preliminary distribution parameters; if the target step is executed a preset number of times, use the preliminary distribution parameters when the preset number of times is reached as the target distribution parameters, or if the preliminary distribution parameters converge, use the preliminary distribution parameters as the target distribution parameters.

3. The method according to claim 2, characterized in that The step of constructing the first log-likelihood function of the preset distribution according to the gene expression level of each of the spatial sequencing sites comprises: Constructing the probability expression data of the preset distribution according to the gene expression amount of each of the spatial sequencing sites; Constructing a state variable of the spatial sequencing site according to the probabilistic expression data; wherein the state variable is used to indicate whether the spatial sequencing site belongs to a cell coverage area or a non-cell coverage area; Constructing a distribution variable of the preset distribution according to the probability expression data; wherein the distribution variable is used to indicate that the type of the preset distribution is a cell coverage area type or a non-cell coverage area type; The first log-likelihood function is constructed according to the state variable and the distribution variable.

4. The method according to claim 2, characterized in that: The conditional probability includes the cell region probability, the first variable probability and the second variable probability, and the original distribution parameters include the first distribution parameter and the second distribution parameter; The step of calculating the conditional probability of the latent variable according to the randomly initialized original distribution parameters includes: Calculate the probability that the spatial sequencing site belongs to the cell coverage area according to the preset conditions and the original distribution parameters to obtain the cell area probability; Perform logarithmic calculation according to the first distribution parameter to obtain the first variable probability; Perform logarithmic differentiation calculation according to the second distribution parameter to obtain the second variable probability.

5. The method according to claim 4, characterized in that The maximizing of the first log-likelihood function according to the conditional probability to obtain preliminary distribution parameters includes: Summing the cell region probabilities to obtain a probability sum; Calculating the product sum of the cell region probability and the second variable probability, and performing ratio calculation on the product sum and the accumulated value of the cell region probability to obtain a first probability ratio; A second probability ratio is calculated according to the cell region probability, the first variable probability, and the second variable probability; An expectation function of the first log-likelihood function is calculated, and the expectation function is maximized according to the probability sum, the first probability ratio, and the second probability ratio to obtain the preliminary distribution parameters.

6. The method according to claim 1, characterized in that The calculating of the target distribution parameters corresponding to the gene expression amount of each of the spatial sequencing sites obeying the preset distribution includes: Constructing a second probability density function of the preset distribution according to the gene expression level of each of the spatial sequencing sites; Obtaining observation data of the spatial sequencing site belonging to the cell coverage area in a preset number of times, and constructing a second likelihood function according to the observation data and the second probability density function; Performing logarithmic processing on the second likelihood function to obtain a logarithmic likelihood function; The log-likelihood function is differentiated to obtain the target distribution parameters.

7. The method according to claim 1, characterized in that Determining a first probability that the spatial sequencing site belongs to a cell coverage area according to the target distribution parameter of the adjacent site includes: Calculate the adjacent message data between the adjacent site and the spatial sequencing site according to the target distribution parameter; A first probability that the spatial sequencing site belongs to a cell coverage area is determined based on the neighboring message data.

8. The method according to claim 7, characterized in that The step of calculating the adjacent message data between the adjacent site and the spatial sequencing site according to the target distribution parameter comprises: Calculating the conditional probability that the spatial sequencing site belongs to the cell coverage area according to the target distribution parameters of the spatial sequencing site; Calculating a first potential energy difference of the adjacent message according to the first state data of the adjacent site, the target distribution parameter of the adjacent site, and a preset potential function; wherein the first state data indicates whether the adjacent site belongs to a cell coverage area or a non-cell coverage area; The adjacent message data is calculated according to the conditional probability and the first potential energy difference.

9. The method according to claim 7, characterized in that: The determining, according to the adjacent message data, a first probability that the spatial sequencing site belongs to a cell coverage area comprises: Summing a plurality of adjacent message data to obtain total message data; Calculating a second potential energy difference of the spatial sequencing site according to the second state data of the spatial sequencing site and the target distribution parameter of the spatial sequencing site; wherein the second state data indicates that the spatial sequencing site belongs to a cell coverage area or a non-cell coverage area; A first probability that the spatial sequencing site belongs to a cell coverage area is calculated based on the total message data and the second potential energy difference.

10. The method according to any one of claims 1 to 9, characterized in that: The step of determining the cell region of the target slice according to the first probability and the position data comprises: constructing a cell mask matrix according to the first probability and the position data, and constructing a spatial count matrix according to the gene expression information; The cell area of ​​the target slice is determined according to the cell mask matrix and the spatial count matrix.

11. The method according to claim 10, characterized in that The step of constructing a cell mask matrix according to the first probability and the position data comprises: Determining a second probability that the spatial sequencing site belongs to a non-cell coverage area according to the target distribution parameters of the adjacent sites; The first probability is compared with the second probability, and the cell mask matrix is ​​constructed according to the comparison result, the preset mask value, and the position data.

12. The method according to claim 10, characterized in that The root determines the cell area of ​​the target slice according to the cell mask matrix and the spatial count matrix, including: Constructing a convolution kernel according to a preset cell diameter, and performing convolution calculation according to the convolution kernel and the spatial counting matrix to obtain a convolution matrix; Performing dot product calculation on the convolution matrix and the cell mask matrix to obtain a cell label matrix; wherein the row and column dimensions of the cell label matrix are the same as those of the spatial count matrix, the cell label matrix includes the position data of the spatial sequencing sites and the cell labels corresponding to the spatial sequencing sites, and the cell labels include cell region labels; The cell region of the target slice is determined according to the position data corresponding to the cell region label.

13. A cell region determination device, characterized in that: The device comprises: An information acquisition unit, used to acquire gene expression information from the spatial transcription data of the target slice; wherein the gene expression information includes position data of multiple spatial sequencing sites and the gene expression level of each of the spatial sequencing sites; A distribution parameter calculation unit, used to calculate the target distribution parameter corresponding to the gene expression amount of each of the spatial sequencing sites when it follows a preset distribution; A probability calculation unit, used to calculate a first probability that each of the spatial sequencing sites belongs to a cell coverage area according to the target distribution parameters of the adjacent sites of each spatial sequencing site; wherein the adjacent sites of the spatial sequencing sites are spatial sequencing sites adjacent to the spatial sequencing sites; A cell region determining unit is used to determine the cell region of the target slice according to the first probability and the position data.

14. An electronic device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the cell region determination method according to any one of claims 1 to 12 is implemented.

15. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the cell region determination method according to any one of claims 1 to 12 is implemented.

16. A computer program product, comprising a computer program, wherein the computer program is read and executed by a processor of a computer device, so that the computer device executes the cell region determination method according to any one of claims 1 to 12.

Images