Multi-slice registration method and device, electronic device and storage medium

Abstract

The present disclosure provides a multi-slice registration method and device, an electronic device and a storage medium. The multi-slice registration method includes: obtaining a first feature set of the source slice based on the spatial transcription data of the source slice, and obtaining a second feature set of the target slice based on the spatial transcription data of the target slice; determining feature difference data according to the first feature set and the second feature set; determining first distribution data and second distribution data; constructing a target transmission plan according to the feature difference data, the distribution difference between the first distribution data and the second distribution data, and registering the source slice and the target slice according to the target transmission plan. The embodiments of the present application can improve the accuracy of 3D registration of unbalanced data.

Description

Technical Field

[0001] The present disclosure relates to the field of bioinformatics, and in particular to a multi-slice registration method and device, an electronic device, and a storage medium. Background Art

[0002] Spatial transcriptomics (ST) is an emerging technology that can measure mRNA expression at thousands of locations in tissue samples while maintaining the original position of cells.

[0003] Currently, a variety of 3D registration software and algorithms have been developed to integrate ST data from different slices in order to reconstruct a high-precision 3D tissue model. However, when the number of points in two slices is different, the methods in the related art will force the slices to match, resulting in biological interpretation errors and low 3D registration accuracy. Therefore, how to provide a multi-slice registration method to improve the accuracy of 3D registration of unbalanced data has become a technical problem that needs to be solved urgently. Summary of the invention

[0004] The main purpose of the embodiments of the present application is to propose a multi-slice registration method and device, an electronic device and a storage medium, aiming to improve the accuracy of 3D registration of unbalanced data.

[0005] To achieve the above object, a first aspect of an embodiment of the present application proposes a multi-slice registration method, the method comprising:

[0006] Acquire a first feature set of the source slice based on the spatial transcription data of the source slice, and acquire a second feature set of the target slice based on the spatial transcription data of the target slice;

[0007] Determining feature difference data according to the first feature set and the second feature set;

[0008] Determine first distribution data and second distribution data, wherein the first distribution data is used to describe the total mass transmitted from the source site in the source slice to the target site in the target slice, and the second distribution data is used to describe the total mass received by the target site;

[0009] A target transmission plan is constructed according to the feature difference data and the distribution difference between the first distribution data and the second distribution data, and the source slice and the target slice are aligned according to the target transmission plan.

[0010] In some embodiments, the first feature set includes first gene expression data and first site coordinate data, and the second feature set includes second gene expression data and second site coordinate data;

[0011] The determining feature difference data according to the first feature set and the second feature set includes:

[0012] determining a first gene expression difference according to the first gene expression data and the second gene expression data;

[0013] Determine a first point distance difference according to the first point coordinate data and the second point coordinate data;

[0014] Determine a first site weight matrix of the source slice and a second site weight matrix of the target slice respectively; the first site weight matrix is ​​used to describe the relative importance of each source site in the source slice; the second site weight matrix is ​​used to describe the relative importance of each target site in the target;

[0015] Determining a first site weight difference according to the first site weight matrix and the second site weight matrix;

[0016] The feature difference data is determined according to the first gene expression difference, the first site distance difference and the first site weight difference.

[0017] In some embodiments, the first feature set includes first gene expression data and first site coordinate data, and the second feature set includes second gene expression data and second site coordinate data;

[0018] The determining feature difference data according to the first feature set and the second feature set includes:

[0019] Sampling the source site to obtain a source sampling site, and sampling the target site to obtain a target sampling site;

[0020] Determine a second gene expression difference according to the source sampling site, the target sampling site, the first gene expression data, and the second gene expression data;

[0021] Determine a second site distance difference according to the source sampling site, the target sampling site, the first site coordinate data, and the second site coordinate data;

[0022] Determine a third site weight matrix of the source slice and a fourth site weight matrix of the target slice; the third site weight matrix is ​​used to describe the relative importance of each source site in the source slice; the fourth site weight matrix is ​​used to describe the relative importance of each target site in the target;

[0023] Determining a second site weight difference according to the third site weight matrix and the fourth site weight matrix;

[0024] The feature difference data is determined according to the second gene expression difference, the second site distance difference and the second site weight difference.

[0025] In some embodiments, before determining feature difference data according to the first feature set and the second feature set, the method further includes performing dimensionality reduction processing on the first gene expression data, including:

[0026] performing normalization processing on the first gene expression data and the second gene expression data respectively;

[0027] Merging the normalized first gene expression data and the second gene expression data according to a common gene dimension to obtain gene merged data;

[0028] Performing dimensionality reduction processing on the gene merged data to obtain total gene dimensionality reduction data;

[0029] Sub-gene dimensionality reduction data of the first gene expression data are extracted from the total gene dimensionality reduction data based on the cell labels of the source slice, and the sub-gene dimensionality reduction data are used as the dimensionality reduction result of the first gene expression data.

[0030] In some embodiments, the first site coordinate data is used to describe the first initial coordinates of each source site of the source slice in the initial space, and the second site coordinate data is used to describe the second initial coordinates of each target site of the target slice in the initial space;

[0031] Before determining feature difference data according to the first feature set and the second feature set, the method further includes performing dimensionality reduction processing on the first site coordinate data, including:

[0032] Determine a site distance matrix according to the first initial coordinates and the second initial coordinates;

[0033] Performing dimensionality reduction processing on the site distance matrix to obtain a first initial low-dimensional coordinate corresponding to each source site of the source slice in the low-dimensional space;

[0034] Determine the nearest neighbor distance between the first initial coordinate of each of the source sites and the first initial low-dimensional coordinate, and determine a reference distance based on a plurality of the nearest neighbor distances;

[0035] The low-dimensional space is scaled according to the reference distance, the first target low-dimensional coordinates of each source site are obtained according to the scaled low-dimensional space, and the dimensionality reduction result of the first site coordinate data is obtained according to the first target low-dimensional coordinates.

[0036] In some embodiments, sampling the source site to obtain the source sampling site includes:

[0037] Determining a segmentation size according to the reference distance and a preset scaling factor;

[0038] Segmenting the low-dimensional space according to the segmentation size to obtain segmentation units;

[0039] The source sampling point is obtained by sampling the source point according to the segmentation unit.

[0040] In some embodiments, constructing a target transmission plan according to the distribution difference between the feature difference data, the first distribution data, and the second distribution data includes:

[0041] constructing an initial transmission plan according to the characteristic difference data and the distribution difference between the first distribution data and the second distribution data;

[0042] Site registration mask data is constructed based on the initial transmission plan, and the target transmission plan is constructed based on the site registration mask data.

[0043] In some embodiments, the initial transmission plan is used to describe the quality of transmission from the source site to the target site;

[0044] The constructing site registration mask data according to the initial transmission plan includes:

[0045] Constructing a cluster set according to the maximum transmission quality in the initial transmission plan, constructing first cluster data of the source slice according to the relationship between the source sampling site and the cluster set, and constructing second cluster data of the target slice according to the relationship between the target sampling site and the cluster set;

[0046] Constructing first neighborhood data of the source slice and second neighborhood data of the target slice according to the cluster set and a preset distance;

[0047] The site registration mask data is constructed according to the first clustering data, the second clustering data, the first neighborhood data, and the second neighborhood data.

[0048] In some embodiments, registering the source slice and the target slice according to the target transmission plan includes:

[0049] Acquire first centroid data of the source site and second centroid data of the target site according to the target transmission plan;

[0050] The registration coordinates are calculated according to the first centroid data and the second centroid data, and the source slice and the target slice are registered according to the registration coordinates.

[0051] To achieve the above-mentioned purpose, a second aspect of an embodiment of the present application provides a multi-slice registration device, the device comprising:

[0052] A feature acquisition unit, configured to acquire a first feature set of a source slice based on the spatial transcription data of the source slice, and acquire a second feature set of a target slice based on the spatial transcription data of the target slice;

[0053] a feature difference determining unit, configured to determine feature difference data according to the first feature set and the second feature set;

[0054] a distribution data determination unit, configured to determine first distribution data and second distribution data, wherein the first distribution data is used to describe the total quality transmitted from a source site in a source slice to a target site in a target slice, and the second distribution data is used to describe the total quality received by the target site;

[0055] The registration unit is used to construct a target transmission plan according to the feature difference data and the distribution difference between the first distribution data and the second distribution data, and to register the source slice and the target slice according to the target transmission plan.

[0056] 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.

[0057] 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.

[0058] 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.

[0059] The multi-slice registration method and device, electronic device and storage medium proposed in the embodiment of the present application determine the feature difference data through the first feature set of the source slice and the second feature set of the target slice, and construct a target transmission plan based on the difference between the feature difference data, the first distribution data and the second distribution data. It can be seen that the target transmission plan of the embodiment of the present application not only takes into account the difference in features between the two slices, but also takes into account the difference between the transport quality of the source slice and the receiving quality of the target slice, that is, the embodiment of the present application allows partial mass transmission. In other words, the embodiment of the present application can handle the registration of spatially unbalanced data. In this way, when the source slice and the target slice are registered according to the target transmission plan, the registration operation is more in line with the actual situation, the situation where forced matching leads to biological interpretation errors is reduced, and the registration accuracy is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] 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.

[0061] Figure 1 is a flow chart of a multi-slice registration method provided in an embodiment of the present application;

[0062] Figures 2 to 3 is a schematic diagram of slice registration provided in an embodiment of the present application;

[0063] Figure 4 is a flow chart of preprocessing the gene expression matrix provided in the embodiment of the present application;

[0064] Figure 5 is a schematic diagram of merging gene expression data provided in an embodiment of the present application;

[0065] Figure 6 is a flow chart of preprocessing the spatial coordinate matrix provided in an embodiment of the present application;

[0066] Figure 7 yes Figure 1 Another embodiment flow chart of step S102 in FIG.

[0067] Figure 8 is a flow chart of a method for sampling a source sampling site provided in an embodiment of the present application;

[0068] Fig. 9 is a flow chart of a method for constructing a target transmission plan provided in an embodiment of the present application;

[0069] Fig.10 yes Fig. 9 Flowchart of step S902 in FIG.

[0070] Fig.11 The embodiment of the present application provides a method for determining unbalanced transmission parameters. Flowchart of the method for obtaining values;

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

[0072] FIG. 13A to FIG. 13I It is a schematic diagram of the effect of adjacent slice registration provided by the application embodiment;

[0073] FIG. 14A to FIG. 14D It is a schematic diagram of the effect of long-distance slice registration provided by the application embodiment;

[0074] FIG. 15A to FIG. 15C is a schematic diagram of a 3D reconstruction effect provided by an embodiment of the application;

[0075] Fig.16 is a schematic diagram of a multi-slice registration device provided in an embodiment of the present application;

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

[0077] 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.

[0078] 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:

[0079] 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 one experiment.

[0080] KL divergence (Kullback-Leibler divergence): also known as relative entropy, is used to measure the difference between two probability distributions. KL divergence is an asymmetric measure of the difference between two probability distributions (such as P and Q), where P is usually considered to be the true distribution and Q is a model distribution or an alternative distribution.

[0081] Highly Variable Genes (HVG): refers to genes whose expression levels vary the most between different cells. These genes have significant expression differences between cells, for example, these genes are highly expressed in some cells and lowly expressed in other cells.

[0082] Principal Component Analysis (PCA): is a data dimensionality reduction technique. PCA transforms variables that may be correlated into a set of linearly uncorrelated variables through orthogonal transformation. The goal of PCA is to reduce the dimensionality of the data while retaining the original data information as much as possible.

[0083] Multi-Dimensional Scaling (MDS): is a statistical method used to transform complex, high-dimensional similarity or distance data into intuitive, low-dimensional visual representations. MDS attempts to find a point matrix in a low-dimensional space so that the distance between each point in the matrix is ​​as close as possible to the similarity or distance measure between the corresponding objects in the original data.

[0084] Singular Value Decomposition (SVD): It is a method of matrix decomposition, also known as matrix factor decomposition, which represents the original matrix as a new structure or the product of two or more matrices with special properties.

[0085] Time complexity: refers to the time required for an algorithm to execute. Time complexity is usually expressed using the parameter O notation, for example: O(1): constant time complexity, indicating that the execution time of the algorithm does not change with the size of the input data. O(n): linear time complexity, indicating that the execution time of the algorithm is proportional to the size of the input data.

[0086] Space complexity: refers to the maximum memory space required during the execution of an algorithm. The representation of space complexity is usually expressed using the parameter O representation, for example: O(1): constant space complexity, indicating that the memory space required by the algorithm during execution does not change with the size of the input data. O(n): linear space complexity, indicating that the memory space required by the algorithm during execution is proportional to the size of the input data.

[0087] In the field of bioinformatics, understanding and describing the detailed distribution and functional activities of cells in tissues is crucial to revealing disease mechanisms and normal physiological processes. Although traditional transcriptomics can reveal the molecular characteristics of cells, it ignores the spatial positioning of cells in tissues. Spatial positioning information is of irreplaceable value for understanding cell-to-cell interactions and signal transduction. The rise of ST technology fills this gap. ST technology can measure mRNA expression at thousands of locations in tissue samples while maintaining the original position of cells. However, ST technology faces multiple challenges when dealing with three-dimensional tissue structures, especially when a complete 3D model needs to be constructed from multiple continuous slices, accurate slice registration becomes the key.

[0088] Currently, a variety of 3D registration software and algorithms have been developed to integrate ST data from different slices in order to reconstruct high-precision 3D tissue models. These tools usually adopt different strategies, such as image feature-based registration, gene expression-based registration, and methods that combine image features and gene expression. For example, image processing toolkits (such as VisuAlign, 3DSlicer) use image registration technology to align images generated by ST slices with the standard brain atlas (CCFv3). The PASTE algorithm is based on the optimal transfer objective function (Fused Gromov-Wasserstein, FGW) and combines expression similarity and spatial similarity for slice registration. The PASTE algorithm is suitable for slice alignment of equal points. The SPACEL algorithm combines deep learning and optimization strategies to register by identifying and comparing similar spatial expression patterns in adjacent slices. Among them, image processing toolkits usually require manual selection of landmark points, which increases the complexity of the operation. In addition, image processing toolkits rely on corresponding histological images, which may limit this method due to cost or equipment availability. The PASTE algorithm forces matching when the number of points in two slices is different, resulting in biological interpretation errors. In addition, the PASTE algorithm is limited by equipment availability and cannot process large-scale data. The accuracy of the SPACEL algorithm depends on accurate spatial region identification and the integrity of the single-cell reference dataset. Incorrect region selection may cause the registration results to deviate from reality.

[0089] In summary, although the above methods have improved the accuracy of 3D registration to a certain extent, these methods generally have the following limitations:

[0090] 1. Most algorithms do not work well when dealing with slices with non-equivalent points, because these algorithms assume that each point of the source slice should completely correspond to a point on the target slice, but this assumption is usually not valid in practical applications, leading to errors in biological interpretation.

[0091] 2. Some tools rely too much on high-quality histological images or the selection of manual landmarks, which increases the complexity and subjectivity of the operation and may also introduce bias due to human factors.

[0092] 3. When the slice interval is large or the tissue structure is complex, the performance of the methods in the relevant technology will be significantly reduced, especially when performing cross-time point analysis, the accuracy and reliability of the registration will face challenges.

[0093] Based on this, the embodiments of the present application provide a multi-slice registration method and device, an electronic device and a storage medium, which can handle the registration problem of spatially unbalanced data and reduce user interference to improve the accuracy and stability of 3D registration.

[0094] The multi-slice registration method provided in the embodiment of the present application relates to the field of bioinformatics. The multi-slice registration method provided in the embodiment of the present application can be applied to a terminal, can be applied to a server side, or can be software running in a terminal or a server side. In some embodiments, the terminal can be a smart phone, a tablet computer, a laptop computer, a desktop computer, etc.; the server side can be configured as an independent physical server, or a server cluster or a distributed system composed of multiple physical servers, or a cloud server that provides basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communications, middleware services, domain name services, security services, CDN, and big data and artificial intelligence platforms; the software can be an application that implements the multi-slice registration method, etc., but is not limited to the above forms.

[0095] The present application can be used in many general or special computer system environments or configurations. For example: personal computers, server computers, handheld or portable devices, tablet devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable consumer electronics, network PCs, minicomputers, mainframe computers, distributed computing environments including any of the above systems or devices, etc. The present application can be described in the general context of computer-executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform specific tasks or implement specific abstract data types. The present application can also be practiced in distributed computing environments, in which tasks are performed by remote processing devices connected through a communication network. In a distributed computing environment, program modules can be located in local and remote computer storage media including storage devices.

[0096] The multi-slice registration method provided in the embodiment of the present application is described below.

[0097] Reference Figure 1 In some embodiments, the multi-slice registration method provided in the embodiments of the present application includes but is not limited to steps S101 to S104.

[0098] Step S101, acquiring a first feature set of a source slice based on spatial transcription data of the source slice, and acquiring a second feature set of a target slice based on spatial transcription data of the target slice;

[0099] Step S102, determining feature difference data according to the first feature set and the second feature set;

[0100] Step S103, determining first distribution data and second distribution data, wherein the first distribution data is used to describe the total quality transmitted from the source site in the source slice to the target site in the target slice, and the second distribution data is used to describe the total quality received by the target site;

[0101] Step S104, constructing a target transmission plan according to the feature difference data and the distribution difference between the first distribution data and the second distribution data, and registering the source slice and the target slice according to the target transmission plan.

[0102] In steps S101 to S104 shown in the embodiment of the present application, feature difference data is determined through the first feature set of the source slice and the second feature set of the target slice, and a target transmission plan is constructed based on the difference between the feature difference data, the first distribution data, and the second distribution data. It can be seen that the target transmission plan of the embodiment of the present application not only takes into account the difference in features between the two slices, but also takes into account the difference between the transport quality of the source slice and the receiving quality of the target slice, that is, the embodiment of the present application allows partial mass transmission. In other words, the embodiment of the present application can handle the registration of spatially unbalanced data. In this way, when the source slice and the target slice are registered according to the target transmission plan, the registration operation is more in line with the actual situation, the situation where forced matching leads to biological interpretation errors is reduced, and the registration accuracy is improved.

[0103] In step S101 of some embodiments, the source slice and the target slice may refer to slices from the same target sample to be registered. The slice may refer to a biological slice obtained in accordance with relevant legal provisions, including plant slices and animal slices. Registration may refer to matching the source site (i.e., source spot) in the spatial transcription data of the source slice with the target site in the spatial transcription data of the target slice to achieve 3D reconstruction based on multi-slice registration. For example, by Figure 2 By registering the multiple slices shown in Figure 3 The registration effect shown ( Figure 3 The black outline in the figure is the real biological position of the slice, and the slice itself is the position after registration). Subsequent 3D reconstruction can be performed based on the registration effect. The spatial transcription data of the source slice can be obtained through ST technology, and the spatial transcription data of the target slice can be obtained. The first feature set can refer to the expression set of the feature elements on the source slice, and the second feature set can refer to the expression set of the feature elements on the target slice. The first feature set can be determined according to the spatial transcription data corresponding to the source slice, and the second feature set can be determined according to the spatial transcription data corresponding to the target slice.

[0104] For example, feature elements may include a gene expression matrix (Indicated by genes and A matrix consisting of sites, where Indicates The gene in expression level at each site), spatial coordinate matrix (in, Indicates The two-dimensional coordinates of the sites, express A real number matrix with 2 rows and 2 columns), a site distance matrix (in, ). Since the placement and orientation of slices in two-dimensional space can be arbitrarily defined, the site distance matrix Eliminate this randomness.

[0105] It is understood that in order to define the relative importance of each site, the feature element can also include a weight matrix , the weight matrix can be used to represent the importance of each site relative to other sites. , this site The corresponding weight It can be determined based on prior knowledge, such as calculated based on specific cell labels or pathological results. After weight normalization, When no prior knowledge is available, it can be assumed that each site has equal weight, such as In this case, the weight of each site can be adjusted according to the subsequent alignment results.

[0106] In step S102 of some embodiments, the feature difference data may refer to data determined based on the expression difference of the same feature element in the first feature set and the second feature set. As the feature element in the above example, the feature difference data may include difference data related to gene expression and difference data related to site distance.

[0107] In step S103 of some embodiments, the first distribution data may refer to a first marginal distribution, which may be used to describe the total mass transmitted from a source site in a source slice to a target site in a target slice. The second distribution data may refer to a second marginal distribution, which may be used to describe the total mass received by the target site in the target slice. It is understood that in the embodiments of the present application, the mass may refer to the probability of a site participating in the transmission.

[0108] In step S104 of some embodiments, a target transmission plan is constructed based on the feature difference data, the distribution difference between the first distribution data and the second distribution data, and a transportation plan can be obtained by minimizing the target transmission plan, and the source slice and the target slice can be registered according to the transportation plan. The transportation plan can be a matrix for describing the transmission probability between any source site and any target site.

[0109] It can be understood that the first feature set can be based on genes and An undirected graph constructed from source slices of source points Similarly, the second feature set can be based on the inclusion of genes and An undirected graph constructed by target slices of target sites Based on this, we can construct the transmission plan shown in the following formula 1 (i.e., the objective function ).

[0110] (Formula 1)

[0111] in, and The meaning is similar to , and The meaning is similar to . The edge constraints shown in Formula 2 below are satisfied.

[0112] (Formula 2)

[0113] It can be understood that in formula 1, the objective function It includes a gene expression similarity term (i.e., the first summation term) and a spatial expression similarity term (i.e., the second summation term). These two terms are connected by parameters Weighted. Among them, the gene expression similarity term can represent the difference in gene expression from each site To each site The cost of moving a unit probability mass. The spatial expression similarity term is used to maintain the spatial distance between sites within the slice. However, in the method shown in Formula 1, under the hard edge constraint, A rigid structure must be maintained, corresponding to the source slice Any source point in All weights must be sent, and for the target slice Any target site in must accept all weights. This constraint applies to the gene expression matrix , site distance matrix This constraint is valid for similar source and target slices, but for source and target slices that are significantly different (such as large morphological differences or loss during sample collection), this constraint will result in only the source slice existing. Source loci corresponding to specific cell types or tissue regions in the target slice are incorrectly mapped to the target slice. Any target site in .

[0114] In some embodiments, step S102 includes but is not limited to the following steps:

[0115] determining a first gene expression difference based on the first gene expression data and the second gene expression data;

[0116] Determine the first point distance difference according to the first point coordinate data and the second point coordinate data;

[0117] Determine the first site weight matrix of the source slice and the second site weight matrix of the target slice respectively; the first site weight matrix is ​​used to describe the relative importance of each source site in the source slice; the second site weight matrix is ​​used to describe the relative importance of each target site in the target;

[0118] Determine the first site weight difference according to the first site weight matrix and the second site weight matrix;

[0119] The feature difference data is determined according to the first gene expression difference, the first site distance difference and the first site weight difference.

[0120] It is understood that in some embodiments, the first gene expression data may be used to and second gene expression data Determine the first gene expression difference. According to the first site coordinate data The site distance matrix can be determined , according to the coordinate data of the second site The site distance matrix can be determined , according to the site distance matrix and site distance matrix The first point distance difference can be determined. According to the first point weight matrix and the second site weight matrix The first site weight difference can be determined. The feature difference data can be determined according to the first gene expression difference, the first site distance difference and the first site weight difference.

[0121] Based on the above characteristic difference data, the embodiment of the present application can construct the target transmission plan (i.e., the target function) shown in the following formula 3: ).

[0122] (Formula 3)

[0123] In formula 3, the definitions of the first two items are the same as those in formula 1, that is, the first item is the first gene expression difference, and the second item is the first site distance difference. In the third item, represents the first distribution data, represents the second distribution data. Different from Formula 1, the target transmission plan shown in Formula 3 measures the difference between the marginal distribution of the current transportation calculation (such as the first distribution data) and the expected distribution (such as the second distribution data) through KL divergence. By minimizing these KL divergence terms, the first distribution data is close to but different from the second distribution data, thereby achieving partial quality transmission. It can be understood that the parameter Used to compare the first distribution data (or second distribution data) with the prior weights (or prior weights ), thereby adjusting the size of the region with significant signal or geometric differences between the source and target slices to reduce false matches. The last term in Equation 3 is used to smooth the calculation process and to speed up the approximate solution of the optimization problem. In practice, the last term is usually combined with a small parameter Combined calculations are performed to reduce the impact on calculation accuracy. The first site weight difference data between the first site weight matrix and the second site weight matrix can be determined according to the KL divergence.

[0124] The embodiment of the present application relaxes the edge conditions so that the transport plan is allowed not to transmit all the mass, that is, each source site in the source slice can transmit less than the mass of the source site itself, and each target site in the target slice can receive more or less than the mass of the target site itself. In other words, under the initial conditions, it can be considered that the probability of each site participating in the transmission is the same, so after normalization, the total mass is 1. In the embodiment of the present application, each site can be allowed to transmit less than its own mass, that is, the probability of the site participating in the transmission is less than or equal to 1. Based on this, the embodiment of the present application realizes more flexible alignment and reduces the problems of incorrect alignment and non-compliance with actual biological conditions in the related art.

[0125] It is understandable that due to the inherent limitations of ST technology in measuring mRNA expression, many sites in the data set lack sufficient gene expression information, which affects the accuracy and completeness of spatial transcription data analysis. Therefore, before step S102, the data can also be preprocessed, that is, the undirected graph corresponding to the source slice can be determined based on the preprocessed data. The undirected graph corresponding to the target slice , in order to reduce the impact of data sparsity on the target transmission plan shown in Formula 3. As follows, the gene expression matrix of the source slice is and the space coordinate matrix It is understandable that the method of preprocessing the target slice is similar to the method of preprocessing the source slice, and this embodiment of the present application will not be described in detail.

[0126] First, the gene expression matrix Refer to Figure 4 In some embodiments, before step S102, the method provided in the embodiment of the present application may further include performing dimensionality reduction processing on the first gene expression data, specifically including but not limited to steps S401 to S404.

[0127] Step S401, normalizing the first gene expression data and the second gene expression data respectively;

[0128] Step S402, merging the normalized first gene expression data and the second gene expression data according to the common gene dimension to obtain gene merged data;

[0129] Step S403, performing dimensionality reduction processing on the gene merged data to obtain total gene dimensionality reduction data;

[0130] Step S404: extract sub-gene dimensionality reduction data of the first gene expression data from the gene dimensionality reduction data based on the cell labels of the source slice, and use the sub-gene dimensionality reduction data as the dimensionality reduction result of the first gene expression data.

[0131] In step S401 of some embodiments, the first gene expression data may refer to the gene expression matrix of the source slice , the second gene expression data may refer to the gene expression matrix of the target slice . Normalization processing can refer to logarithmic normalization processing and scaling processing. Logarithmic normalization processing can be used to compress the data range to a smaller interval. When the original data has maximum and minimum values, or the distribution is uneven, the use of logarithmic normalization can reduce the impact of these extreme values ​​on subsequent analysis, making the data more concentrated in the central area. Scaling processing can refer to adjusting the data in the matrix to a specific range to facilitate comparison and processing between different features. The first gene expression data and the second gene expression data are normalized separately.

[0132] In step S402 of some embodiments, the data structure of the normalized first gene expression data and the second gene expression data can be a DataFrame matrix (DataFrame matrix is ​​one of the data structures in pandas, similar to a matrix or table format) of pandas (pandas is a data analysis and operation library, pandas provides a variety of data structures to process and analyze data). The row labels of the normalized first gene expression data can be the cell labels or cell names of the source slice, and the column labels are the genes shared by the source slice and the target slice. Similarly, the row labels of the normalized second gene expression data can be the cell labels or cell names of the target slice, and the column labels are the genes shared by the target slice and the source slice. The normalized first gene expression data and the second gene expression data are merged according to the shared gene dimension to obtain gene merged data. It can be understood that the number of columns in the gene merged data remains unchanged, but the number of rows increases. For example, as Figure 5 As shown, the normalized first gene expression data 501 and the normalized second gene expression data 502 are merged based on the common genes (ie, gene 1 and gene 2) to obtain gene merged data 503.

[0133] In step S403 of some embodiments, the gene merge data is subjected to dimensionality reduction data to obtain total gene dimensionality reduction data. It is understandable that the dimensionality reduction processing method can be adaptively set according to actual needs, and this embodiment of the present application is not specifically limited to this. For example, the highly variable genes in the gene merge data can be determined by means of Seurat's FindVariableFeatures function (Seurat is a single-cell data analysis tool), Monocle (Monocle is an R package for single-cell transcriptome analysis), etc. The determined highly variable genes are sorted, and the gene merge data is subjected to dimensionality reduction processing based on the highly variable genes ranked in the top N positions to obtain total gene dimensionality reduction data. It is understandable that the specific value of N can be adaptively set according to actual needs, and this embodiment of the present application is not specifically limited to this. For example, the top 3000 highly variable genes can be selected to perform PCA dimensionality reduction to 50 dimensions on the gene merge data to obtain total gene dimensionality reduction data.

[0134] In step S404 of some embodiments, corresponding data (ie, sub-gene dimensionality reduction data) is extracted from the gene dimensionality reduction data based on the cell labels of the source slice, and the extracted data is used as the dimensionality reduction result of the first gene expression data.

[0135] The benefit of step S401 to step S404 is that in the spatial transcriptomics technology based on sequencing, the gene expression matrix is ​​usually a high-dimensional sparse matrix, and the gene expression matrix may contain tens of thousands of genes, while only a few genes are expressed in each cell. By pre-processing the gene expression matrix to reduce the dimension, the complexity of the subsequent calculation transmission plan can be reduced. It is understandable that in some embodiments, other methods (such as singular value decomposition, etc.) can also be used to pre-process the gene expression matrix, which is not specifically limited in the embodiments of the present application.

[0136] Secondly, the space coordinate matrix Refer to Figure 6 In some embodiments, before step S102, the method provided in the embodiment of the present application may further include performing dimensionality reduction processing on the first point coordinate data, specifically including but not limited to steps S601 to S604.

[0137] Step S601, determining a site distance matrix according to the first initial coordinates and the second initial coordinates;

[0138] Step S602, performing dimensionality reduction processing on the site distance matrix to obtain a first initial low-dimensional coordinate corresponding to each source site of the source slice in the low-dimensional space;

[0139] Step S603, determining the nearest neighbor distance between the first initial coordinate of each source point and the first initial low-dimensional coordinate, and determining a reference distance according to a plurality of nearest neighbor distances;

[0140] Step S604, scaling the low-dimensional space according to the reference distance, obtaining the first target low-dimensional coordinates of each source point according to the scaled low-dimensional space, and obtaining the dimensionality reduction result of the first point coordinate data according to the first target low-dimensional coordinates.

[0141] In step S601 of some embodiments, the first point coordinate data may refer to the spatial coordinate matrix of the source slice , the second site coordinate data can refer to the spatial coordinate matrix of the target slice . The first initial coordinate may refer to an element in the first position coordinate data, and the first initial coordinate may refer to the coordinate of the source position of the source slice in the initial space (i.e., the original space, which refers to the natural state of the data before any dimensionality reduction processing). Similarly, the second initial coordinate may refer to the coordinate of the target position of the target slice in the initial space. The position distance matrix may refer to a matrix used to describe the spatial distance between the source position and the target position. For example, the position distance matrix may be a distance matrix as shown in the following formula 4: .

[0142] (Formula 4)

[0143] in, Indicates The first initial coordinates of the source points, Indicates The second initial coordinates of the target site.

[0144] In step S602 of some embodiments, the dimensionality reduction process may refer to multidimensional scaling analysis. By performing multidimensional scaling analysis on the site distance matrix, the source slices can be obtained in the low-dimensional space. The first initial low-dimensional coordinates of the source sites and the first initial low-dimensional coordinates of the target slice in the low-dimensional space The second initial low-dimensional coordinates of the target site. It can be understood that the distance structure between the first initial low-dimensional coordinates and the second initial low-dimensional coordinates in the low-dimensional space is similar to the distance structure between the first initial coordinates and the second initial coordinates in the initial space.

[0145] Specifically, multidimensional scaling analysis may be performed according to the following formula 5 and formula 6.

[0146] (Formula 5)

[0147] (Formula 6)

[0148] in, represents the centering matrix. Represents the dual center distance matrix, the matrix Used to remove overall bias in the data, focusing on the relative distance structure between points. represents the matrix used to calculate the mean, represents the identity matrix, According to the distance matrix Calculated.

[0149] Understandably, due to the computational complexity of multidimensional scaling analysis The higher the value, the higher the computational cost will be when processing large-scale data (e.g., more than 1,000 sites). Therefore, singular value decomposition can be used to further process the low-dimensional coordinates obtained based on multidimensional scaling analysis to speed up the calculation process and reduce the computational cost. For example, singular value decomposition can be performed according to the following formulas 7 and 8.

[0150] (Formula 7)

[0151] (Formula 8)

[0152] in, is the covariance matrix that captures the variance structure of the data in the doubly centered space. represents the feature vector, represents the eigenvalue matrix. is a matrix The characteristic vector of is the corresponding eigenvalue.

[0153] In step S603 of some embodiments, for any source point, a distance calculation is performed based on the first initial coordinate of the source point and any first initial low-dimensional coordinate of the low-dimensional space to determine the nearest neighbor distance. The median of the nearest neighbor distances corresponding to multiple source points is determined to obtain a reference distance.

[0154] In step S604 of some embodiments, the low-dimensional space is scaled according to the reference distance so that the coordinates in the low-dimensional space have the same scale as the coordinates in the initial space. The coordinates of the source site in the low-dimensional space after the scaling process are determined to obtain the first target low-dimensional coordinates. The spatial coordinate matrix of the source slice can be constructed based on the first target low-dimensional coordinates of the multiple source sites, thereby achieving dimensionality reduction of the first point coordinate data. For example, the low-dimensional space can be scaled according to the following formula 9 and formula 10.

[0155] (Formula 9)

[0156] (Formula 10)

[0157] in, It means that the eigenvectors corresponding to the first two largest eigenvalues ​​are selected as the low-dimensional embedding of the distance matrix, that is, the representation of the distance matrix in the low-dimensional space. is the final low-dimensional embedding coordinate (i.e., the dimensionality reduction result of the first point coordinate data).

[0158] It can be seen that in other embodiments, the original data of each feature in Formula 3 can be replaced with the data after the above preprocessing, so as to reduce the problem of data sparsity affecting the accuracy of the target transmission plan.

[0159] From the above example, we can see that through the source slice ( is the original data or the data after dimension reduction), the target slice ( The target transmission plan shown in Formula 3 is minimized by using the original data or the data after dimension reduction, so that the transportation plan of the source slice and the target slice can be obtained. However, the computational complexity of this method is , when the number of source sites or target sites is too large, the computational cost will increase. For example, when the method of the embodiment of the present application is used for calculation on the GPU, a large amount of video memory will be occupied. Due to the GPU memory limitation, it is impossible to process larger-scale data on the GPU. When calculating on the CPU, although there is no storage limitation, the calculation speed will be affected by the amount of data. Based on this, the embodiment of the present application provides a method for calculating a transportation plan from coarse to fine. First, a rough transportation plan of the sampled part is calculated on the GPU to determine which sites should be considered in the fine mapping based on the rough transportation plan. Then, a sparse mask is constructed based on the rough transportation plan. Finally, the precise transportation plan of the masked part is calculated on the CPU. As follows, this method for calculating a transportation plan from coarse to fine is described.

[0160] Reference Figure 7 In some embodiments, step S102 may include but is not limited to steps S701 to S706.

[0161] Step S701, sampling the source site to obtain the source sampling site, and sampling the target site to obtain the target sampling site;

[0162] Step S702, determining a second gene expression difference according to the source sampling site, the target sampling site, the first gene expression data, and the second gene expression data;

[0163] Step S703, determining the distance difference of the second site according to the source sampling site, the target sampling site, the first site coordinate data, and the second site coordinate data;

[0164] Step S704, determining a third site weight matrix of the source slice and a fourth site weight matrix of the target slice; the third site weight matrix is ​​used to describe the relative importance of each source site in the source slice; the fourth site weight matrix is ​​used to describe the relative importance of each target site in the target;

[0165] Step S705, determining the second site weight difference according to the third site weight matrix and the fourth site weight matrix;

[0166] Step S706, determining feature difference data according to the second gene expression difference, the second site distance difference and the second site weight difference.

[0167] In step S701 of some embodiments, when dimensionality reduction processing is not performed, the source sampling site may refer to a site obtained by sampling the source site of the source slice in the initial space. The target sampling site may refer to a site obtained by sampling the target site of the target slice in the initial space.

[0168] In the case of dimensionality reduction, the source sampling site may refer to a site obtained by sampling the source site of the source slice in the low-dimensional space. The target sampling site may refer to a site obtained by sampling the target site of the target slice in the low-dimensional space.

[0169] It is understandable that the sampling method may be random sampling or other methods, which is not specifically limited in the embodiments of the present application.

[0170] Reference Figure 8 In some embodiments, when performing dimensionality reduction processing, the method of "sampling the source site to obtain the source sampling site" in step S701 may include but is not limited to steps S801 to S803.

[0171] Step S801, determining a segmentation size according to a reference distance and a preset scaling factor;

[0172] Step S802, dividing the low-dimensional space according to the segmentation size to obtain segmentation units;

[0173] Step S803: sampling the source points according to the segmentation unit to obtain source sampling points.

[0174] In step S801 of some embodiments, the segmentation size may refer to the size used to segment the low-dimensional space, such as the segmentation length. The segmentation size may be obtained by multiplying the reference distance and a preset scaling factor. The specific value of the scaling factor may be adaptively set according to actual needs, and the embodiments of the present application do not specifically limit this.

[0175] In step S802 of some embodiments, the low-dimensional space is gridded according to the segmentation size to divide the low-dimensional space into a plurality of square grids (i.e., segmentation units) of the same size. It is understandable that the shape of the segmentation unit can also be adaptively set according to actual needs, and this embodiment of the present application does not specifically limit this.

[0176] In step S803 of some embodiments, a source point is extracted from each segmentation unit as a source sampling point. For example, for each segmentation unit, the source point closest to the center of the segmentation unit can be used as the source sampling point of the segmentation unit. Alternatively, each segmentation unit is sampled in a random sampling manner to obtain a source sampling point corresponding to each segmentation unit.

[0177] It can be understood that the method for determining the target sampling site is similar to the method for determining the source sampling site, and this embodiment of the present application will not be described in detail.

[0178] In step S702 of some embodiments, without performing dimensionality reduction processing, the gene expression of the source sampling site is determined from the first gene expression data, and the gene expression of the target sampling site is determined from the second gene expression data. The second gene expression difference can be determined based on the gene expression of the source sampling site and the gene expression of the target sampling site.

[0179] In the case of dimensionality reduction, the gene expression of the source sampling site is determined from the dimensionality reduction result of the first gene expression data, and the gene expression of the target sampling site is determined from the dimensionality reduction result of the second gene expression data. The second gene expression difference can be determined based on the gene expression of the source sampling site and the gene expression of the target sampling site.

[0180] In step S703 of some embodiments, without dimensionality reduction processing, the distance matrix corresponding to the source slice can be obtained according to the first site coordinate data, and the distance matrix corresponding to the target slice can be obtained according to the second site coordinate data. The second site distance difference is determined according to the expression of the source sampling site in the source slice distance matrix and the expression of the target sampling site in the target slice distance matrix.

[0181] In the case of dimensionality reduction, the low-dimensional embedding of the distance matrix corresponding to the source slice can be obtained according to the dimensionality reduction result of the first site coordinate data, and the low-dimensional embedding of the distance matrix corresponding to the target slice can be obtained according to the dimensionality reduction result of the second site coordinate data. The second site distance difference is determined according to the expression of the source sampling site in the low-dimensional embedding of the source slice distance matrix and the expression of the target sampling site in the low-dimensional embedding of the target slice distance matrix.

[0182] In steps S704 to S705 of some embodiments, the third site weight matrix can be used to describe the importance of each source acquisition site in the source slice compared to other source acquisition sites. The fourth site weight matrix can be used to describe the importance of each target acquisition site in the target slice compared to other target acquisition sites. The second site weight difference can be determined based on the third site weight matrix and the fourth site weight matrix.

[0183] In step S706 of some embodiments, the second gene expression difference, the second site distance difference and the second site weight difference may be used as feature difference data.

[0184] Reference Fig. 9 In some embodiments, the method of "constructing a target transmission plan based on feature difference data, first distribution data, and second distribution data" in step S104 may include but is not limited to steps S901 to S902.

[0185] Step S901, constructing an initial transmission plan according to the characteristic difference data, the distribution difference between the first distribution data and the second distribution data;

[0186] Step S902: constructing site registration mask data according to the initial transmission plan, and constructing a target transmission plan according to the site registration mask data.

[0187] In step S901 of some embodiments, an initial transmission plan shown in the following formula 11 may be constructed based on the feature difference data, the first distribution data and the second distribution data. It can be understood that a coarse transmission plan may be obtained based on the minimized initial transmission plan.

[0188] (Formula 11)

[0189] Among them, the first item in formula 11 can be the second gene expression difference in the feature difference data, the second item can be the second site distance expression difference in the feature difference data, the third item measures the difference between the first distribution data and the second distribution data through KL divergence, and the fourth item is used to speed up the calculation convergence.

[0190] By comparison, it can be seen that the initial transmission plan described by Formula 11 is similar to Formula 3, the difference being that Formula 11 is constructed based on the data corresponding to the source sampling site and the target sampling site.

[0191] In step S902 of some embodiments, the site registration mask data may refer to data used to describe which sites (including source sites and target sites) should be considered in the fine mapping. According to the site registration mask data, accurate transportation calculation is performed on the sites to be considered to obtain a target transportation plan.

[0192] Reference Fig.10 In some embodiments, "constructing site registration mask data according to the initial transmission plan" in step S902 may include but is not limited to steps S1001 to S1003.

[0193] Step S1001, constructing a cluster set according to the maximum transmission quality in the initial transmission plan, constructing first cluster data of the source slice according to the relationship between the source sampling site and the cluster set, and constructing second cluster data of the target slice according to the relationship between the target sampling site and the cluster set;

[0194] Step S1002, constructing first neighborhood data of the source slice and second neighborhood data of the target slice according to the cluster set and the preset distance;

[0195] Step S1003: constructing site registration mask data according to the first cluster data, the second cluster data, the first neighborhood data and the second neighborhood data.

[0196] In step S1001 of some embodiments, the rough transmission plan may be in the form of a matrix, and the dimension of the rough transmission plan is the same as the sampling result. In the matrix corresponding to the rough transmission plan, the row labels represent the source sampling sites, the column labels represent the target sampling sites, and the elements in the matrix are as follows: Indicates The source sampling point is The cluster set can be used to cluster the transmission quality with the largest value in each row of the matrix, and the transmission quality with the largest value in each column. For example, assuming that the matrix size corresponding to the rough transmission plan is , then after taking the maximum value of each row and each column respectively, we can get The transmission quality with the largest value. The first cluster data can be used to describe the relationship between each source sampling site and the cluster set. For example, for any source sampling site and any cluster category The clustering data shown in the following formula 12 can be constructed.

[0197] (Formula 12)

[0198] Clustering Data For any source sampling point , when the source sampling point When the point belongs to the cluster set (that is, the source sampling site The transmission quality to a target sampling site belongs to the cluster category ), in the first cluster data The source sampling point is defined as 1, otherwise it is defined as 0. Perform the above operation for each source sampling site to obtain the complete first cluster data Similarly, the second clustering data can be constructed based on the relationship between the target sampling site and the cluster set. .

[0199] In step S1002 of some embodiments, the clustering category (i.e. The first neighborhood data can be used to indicate whether other source sampling sites are within the radius. For example, for any source sampling site (i.e., the source sampling site as the center) and any source sampling site (i.e., other source sampling sites except the center) can be constructed to obtain the neighborhood data as shown in the following formula 13.

[0200] (Formula 13)

[0201] Neighborhood data For any source sampling point , when the source sampling point Source sampling site At the source sampling point Center, distance When the radius is within the range, the first field data The source sampling point is defined as 1, otherwise it is defined as 0. Perform the above operation on the source sampling point corresponding to each center to obtain the complete first field data Similarly, the second clustering data can be constructed .

[0202] In step S1003 of some embodiments, site registration mask data mask shown in the following formula 14 can be constructed based on the first cluster data, the second cluster data, the first neighborhood data and the second neighborhood data.

[0203] (Formula 14)

[0204] It can be understood that in the site registration mask data, the sites defined as 1 represent the sites that need to participate in subsequent transportation, that is, the sites that need to be considered in the fine mapping.

[0205] It can be understood that, according to the site registration mask data mask, an initial transmission plan as shown in the following formula 15 can be obtained.

[0206] (Formula 15)

[0207] According to the initial transmission plan, a target transmission plan as shown in the following formula 16 can be obtained, and a refined transmission plan can be obtained by minimizing the target transmission plan.

[0208] (Formula 16)

[0209] It can be understood that in formula 16, the unbalanced transmission parameter and the fine-grained transmission plan There is a corresponding relationship. When the unbalanced transmission parameter When the value is too large, the entire mass will be forced to be transmitted, resulting in erroneous transmission. If the value is too small, the transmission plan will fail. Therefore, determine the balanced transmission parameters The appropriate value of is crucial. Fig.11 In the embodiment of the present application, the unbalanced transmission parameter can be determined based on the dichotomy method. Specifically, the output result when the transmission plan fails is defined as , to calculate a more accurate refined transmission plan. Get the total quality and , unbalanced transmission parameters The initial value range of and error data . According to the initial value range Calculate unbalanced transmission parameters The initial value of , determine the sum of the fine transmission plan under the initial value state Is it greater than When the detailed transmission plan and Greater than When the upper limit of the initial value range is updated to the unbalanced transmission parameter The initial value of ), in order to narrow the value range. Less than or equal to When the lower limit of the initial value range is updated to the unbalanced transmission parameter The initial value of ). Determine whether the difference between the two end points in the reduced value range is less than the error data (Right now ), when judged as "yes", the unbalanced transmission parameters are output The initial value of and the corresponding fine transmission plan , that is, the initial value is taken as the unbalanced transmission parameter When the judgment is "no", the initial value range can be Update to redefine unbalanced transmission parameters based on new value ranges It is understandable that the initial value range is and error data It can be adaptively set according to actual needs, and this embodiment of the present application does not specifically limit this. The specific value of can also be adaptively set according to actual operation needs.

[0210] Reference Fig.12 As shown, the method of "aligning the source slice and the target slice according to the target transmission plan" in step S104 may include but is not limited to steps S1201 to S1202.

[0211] Step S1201, acquiring first centroid data of a source site and second centroid data of a target site according to a target transmission plan;

[0212] Step S1202: Calculate the registration coordinates according to the first centroid data and the second centroid data, and register the source slice and the target slice according to the registration coordinates.

[0213] In step S1201 of some embodiments, the source slice space coordinate matrix is ​​determined according to the refined transmission plan corresponding to the target transmission plan The centroid of the target slice (i.e., the first centroid data) and the target slice space coordinate matrix The centroid of (i.e., the second centroid data). Specifically, the first centroid data The second centroid data can be calculated according to the following formula 17: It can be calculated according to the following formula 18.

[0214] (Formula 17)

[0215] (Formula 18)

[0216] In step S1202 of some embodiments, the registration coordinates may refer to the original coordinates (i.e., the source slice space coordinate matrix and the target slice space coordinate matrix ) and the refined transmission plan, fixed source slice space coordinate matrix , the target slice space coordinate matrix After performing the translation and rotation without scaling, the new spatial coordinate matrix of the target slice is obtained. It can be seen that determining the registration coordinates is a rotation and translation problem that minimizes the distance between matching sites. This problem can be viewed as solving an orthogonal Procrustes problem. When viewed as solving an orthogonal Procrustes problem, the input coordinates need to have the same length, but usually the lengths of the input coordinates are not the same. To solve this problem, the following formula 19 can be used based on the spatial coordinate matrix Generate coordinates , similarly based on the space coordinate matrix Generate coordinates .

[0217] (Formula 19)

[0218] In formula 19, Represents a refined transmission plan. The refined transmission plan can be in the form of a matrix, where the elements Indicates the source slice source site to target slice The quality of transmission to each target site.

[0219] Referring to the following formula 20, the spatial coordinate matrix is ​​generated according to the sum of the elements in each column of the refined transmission plan: The coordinates of each point in the image are scaled according to the received mass.

[0220] (Formula 20)

[0221] coordinate With coordinates The shapes are the same, both , specifically satisfying the following formula 21.

[0222] (Formula 21)

[0223] Referring to Formula 22, the sampling singular value decomposition method is used to determine the covariance matrix H, where the coordinates and coordinates The normalization of can refer to Formula 22 to Formula 26.

[0224] (Formula 22)

[0225] (Formula 23)

[0226] (Formula 24)

[0227] (Formula 25)

[0228] (Formula 26)

[0229] In Formula 22, are the singular values ​​of the covariance matrix H, and represents an orthogonal matrix. Thus, according to the orthogonal Procrustes problem, we can get , proportionality coefficient In the registration task, there is no need to scale the target slice, so the registration coordinates can be obtained according to the following formula 27 .

[0230] (Formula 27)

[0231] The method provided in the embodiment of the present application performs dimensionality reduction processing on the data by means of PCA, MDS, etc., to ensure the consistency of local distances under rigid body transformation, and to use the gene expression matrix as a feature. In the slice alignment process, the embodiment of the present application provides a coarse-to-fine method. First, a subset of the source slice and the target slice (i.e., the sampling site) is preliminarily calculated on the GPU to generate an initial alignment transmission plan. Then, the transmission plan is refined on the CPU, and the transmission plan is optimized by selecting sites for precise mapping within the fine site registration mask data.

[0232] In a specific embodiment, the embodiment of the present application is verified on multiple data sets such as Barseq and Slideseq. On these data sets, the embodiment of the present application demonstrates advantages in 3D reconstruction and cross-time analysis. Specifically:

[0233] 1. The registration deviation of adjacent slices is extremely small.

[0234] Reference FIG. 13A to FIG. 13I , when processing close, continuous slices, the embodiments of the present application demonstrate excellent registration capabilities. Fig.13A and Fig. 13B Display data examples, N represents the number of slices, and Distance represents the distance between slices. Fig. 13C Indicates that slices are randomly translated and rotated to simulate the actual situation. Fig.13D These are example results of different registration methods (including the FUGW method (i.e., the method of the embodiment of the present application), the PASTE2 method, the SPACEL method, the MOSCOT_Affine method, the MOSCOT_Warp method, and the PASTE method) on the Barseq dataset. Fig.13E It is a box plot of MAE (Mean Absolute Error) on the Barseq dataset. Fig.13F It is a box plot of MJSD (Mean Jensen-Shannon Divergence) on the Barseq dataset. The lower the value, the higher the accuracy. By comparison, it can be seen that the value of the method in the embodiment of the present application (i.e., the FUGW method) is lower than that of other methods, and the Wilcoxon signed rank test is statistically significant. Figure 13G are example results of different registration methods on the Slideseq dataset. Fig.13H It is a box plot of MAE on the Slideseq dataset. Fig.13IIt is a box plot of the MJSD value on the data set Slideseq. By comparison, it can be seen that the value of the method in the embodiment of the present application is lower than that of other methods, and the Wilcoxon signed rank test is statistically significant. By aligning a series of closely arranged slices, the embodiment of the present application can achieve low MAE and MJSD values, which means that the coordinate system after alignment is highly consistent with the actual situation. For example, in the alignment of the 6th slice Slice06 of the data set Barseq as the source slice (as a reference set), the 7th slice Slice07 of the target slice (for translation and rotation transformation) is performed. According to the MJSD distribution diagram, compared with the standard reference coordinates, the deviation of the alignment result of the embodiment of the present application is the smallest, and in some cases, the deviation is almost zero, which is significantly better than other methods (such as PASTE, MOSCOT_Warp and SPACEL).

[0235] 2. It also has a good configuration effect for long-distance slicing.

[0236] Reference FIG. 14A to FIG. 14D For slices that are far apart, the embodiment of the present application also has a better configuration effect. Fig.14A In the data example shown, Distance represents the slice distance. Fig. 14B are example results of different registration methods on the Barseq dataset. Fig. 14C is a box plot of MAE on the Barseq dataset, Fig.14D It is a box plot of the MJSD value on the Barseq dataset. By comparison, it can be seen that the value of the method in the embodiment of the present application is lower than that of other methods, and the Wilcoxon signed rank test is statistically significant. When the slice interval increases to about 900 microns (such as Fig.14A As shown in the figure), it will increase the difficulty of registration because there may be significant anatomical differences between slices. Under such conditions, the embodiment of the present application can still achieve lower average MAE and MJSD values. For example, referring to Fig. 14B As shown in the figure, in the registration of Slice06 as the source slice and Slice11 as the target slice, the performance of the embodiment of the present application is better than that of PASTE2, PASTE, MOSCOT_Affine (abbreviated as Affine) and MOSCOT_Warp (abbreviated as Warp). In particular, compared with PASTE and MOSCOT_Warp, PASTE and MOSCOT_Warp experienced flipping and scaling during registration, resulting in completely wrong results, while the registration results of the embodiment of the present application are more accurate and consistent. It can be understood that in Fig. 14B In the figure, the black outline is the true biological position of the target slice of Slice11. Fig. 14BThe five slice positions shown represent the registration results obtained using different registration methods. The closer the slice position corresponds to the black outline position, the better the registration effect.

[0237] 3. Better 3D reconstruction effect.

[0238] Reference FIG. 15A to FIG. 15C , on the Barseq dataset, the embodiment of the present application achieves better 3D reconstruction results. Fig.15A is the KS test cumulative error curve of 3D reconstruction of different registration methods on the Barseq dataset, Fig. 15B It is the cumulative error sum curve of 3D reconstruction of different registration methods on the Barseq dataset. The lower the value, the higher the accuracy. Fig. 15C It is a comparison of the 3D reconstruction results of different registration methods with Groundtruth (i.e. the real situation). Compared with SPACEL and PASTE2, the cumulative error curve and the error sum curve of the mean absolute error (MAE) of the embodiment of the present application are both statistically significant, with specific p values ​​of 0.027, 5.604×10^-10 (cumulative error curve) and 0.024, 1.868×10^-5 (error sum curve). This shows that the 3D reconstruction results of the embodiment of the present application on the Barseq dataset are better than the other two methods.

[0239] Reference Fig.16 The embodiment of the present application further provides a multi-slice registration device, the device comprising:

[0240] The feature acquisition unit 1601 is used to acquire a first feature set of the source slice based on the spatial transcription data of the source slice, and acquire a second feature set of the target slice based on the spatial transcription data of the target slice;

[0241] A feature difference determining unit 1602, configured to determine feature difference data according to the first feature set and the second feature set;

[0242] The distribution data determining unit 1603 is used to determine first distribution data and second distribution data, wherein the first distribution data is used to describe the total quality transmitted from the source site in the source slice to the target site in the target slice, and the second distribution data is used to describe the total quality received by the target site;

[0243] The registration unit 1604 is used to construct a target transmission plan according to the feature difference data and the distribution difference between the first distribution data and the second distribution data, and to register the source slice and the target slice according to the target transmission plan.

[0244] It can be seen that the contents of the above-mentioned multi-slice registration method embodiment are all applicable to the embodiments of the present multi-slice registration device. The functions specifically implemented by the embodiments of the present multi-slice registration device are the same as those of the above-mentioned multi-slice registration method embodiment, and the beneficial effects achieved are also the same as those achieved by the above-mentioned multi-slice registration method embodiment.

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

[0246] The processor 1701 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;

[0247] The memory 1702 may 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 1702 may store an operating system and other application programs. When the technical solution provided in the embodiment of this specification is implemented by software or firmware, the relevant program code is stored in the memory 1702, and the processor 1701 calls and executes the multi-slice registration method of the embodiment of this application;

[0248] Input / output interface 1703, used to implement information input and output;

[0249] Communication interface 1704, used to realize communication interaction between the device and other devices, which can be realized through wired mode (such as USB, network cable, etc.) or wireless mode (such as mobile network, WIFI, Bluetooth, etc.);

[0250] A bus 1705 that transmits information between various components of the device (e.g., the processor 1701, the memory 1702, the input / output interface 1703, and the communication interface 1704);

[0251] The processor 1701 , the memory 1702 , the input / output interface 1703 and the communication interface 1704 are connected to each other in communication within the device via the bus 1705 .

[0252] 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 multi-slice registration method.

[0253] An embodiment of the present application further provides a computer-readable storage medium, which stores a computer program. When the computer program is executed by a processor, the multi-slice registration method is implemented.

[0254] The memory, as a non-transient computer-readable storage medium, can be used to store non-transient software programs and non-transient computer executable programs. In addition, the memory may include a high-speed random access memory, and may also include a non-transient memory, such as at least one disk storage device, a flash memory device, or other non-transient solid-state storage device. In some embodiments, the memory may optionally include a memory remotely disposed relative to the processor, and these remote memories may be connected to the processor via a network. Examples of the above-mentioned network include, but are not limited to, the Internet, an intranet, a local area network, a mobile communication network, and combinations thereof.

[0255] The embodiments described in the embodiments of the present application are intended to more clearly illustrate the technical solutions of the embodiments of the present application and do not constitute a limitation on the technical solutions provided in the embodiments of the present application. Those skilled in the art will appreciate that with the evolution of technology and the emergence of new application scenarios, the technical solutions provided in the embodiments of the present application are also applicable to similar technical problems.

[0256] Those skilled in the art will appreciate that the technical solutions shown in the figures do not constitute a limitation on the embodiments of the present application, and may include more or fewer steps than shown in the figures, or a combination of certain steps, or different steps.

[0257] The device embodiments described above are merely illustrative, and the units described as separate components may or may not be physically separated, that is, they may be located in one place or distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0258] Those skilled in the art will appreciate that all or some of the steps in the methods disclosed above, and the functional modules / units in the systems and devices may be implemented as software, firmware, hardware, or a suitable combination thereof.

[0259] The terms "first", "second", "third", "fourth", etc. (if any) in the specification of the present application 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 application described herein can be implemented in an order other than those illustrated or described herein. In addition, the terms "including" and "having" and any of their variations are intended to cover non-exclusive inclusions, for example, a process, method, system, product or device comprising a series of steps or units is not necessarily limited to those steps or units clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0260] It should be understood that in the present application, "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.

[0261] In the several embodiments provided in the present application, it should be understood that the disclosed devices and methods can be implemented in other ways. For example, the device embodiments described above are only schematic. For example, the division of the above units is only a logical function division. There may be other division methods in actual implementation, such as 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 through some interfaces, indirect coupling or communication connection of devices or units, which can be electrical, mechanical or other forms.

[0262] The units described above 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.

[0263] In addition, each functional unit in each embodiment of the present application 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.

[0264] 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 application, 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 multiple instructions to enable a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods of various embodiments of the present application. 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 programs.

[0265] The preferred embodiments of the present invention are described above with reference to the accompanying drawings, but the scope of the rights of the present invention is not limited thereto. Any modification, equivalent substitution and improvement made by a person skilled in the art without departing from the scope and essence of the present invention should be within the scope of the rights of the present invention.

Claims

1. A multi-slice registration method, characterized in that: The method comprises: Acquire a first feature set of the source slice based on the spatial transcription data of the source slice, and acquire a second feature set of the target slice based on the spatial transcription data of the target slice; Determining feature difference data according to the first feature set and the second feature set; Determine first distribution data and second distribution data, wherein the first distribution data is used to describe the total mass transmitted from the source site in the source slice to the target site in the target slice, and the second distribution data is used to describe the total mass received by the target site; constructing a target transmission plan according to the feature difference data and the distribution difference between the first distribution data and the second distribution data, and registering the source slice and the target slice according to the target transmission plan; The step of constructing a target transmission plan according to the characteristic difference data and the distribution difference between the first distribution data and the second distribution data includes: constructing an initial transmission plan according to the characteristic difference data and the distribution difference between the first distribution data and the second distribution data; Constructing site registration mask data according to the initial transmission plan, and constructing the target transmission plan according to the site registration mask data; wherein the site registration mask data is used to describe the sites that need to be considered in the fine mapping; The initial transmission plan is used to describe the quality of transmission from the source site to the target site; The constructing site registration mask data according to the initial transmission plan includes: Constructing a cluster set according to the maximum transmission quality in the initial transmission plan, constructing first cluster data of the source slice according to the relationship between the source sampling site and the cluster set, and constructing second cluster data of the target slice according to the relationship between the target sampling site and the cluster set; Constructing first neighborhood data of the source slice and second neighborhood data of the target slice according to the cluster set and a preset distance; The site registration mask data is constructed according to the first clustering data, the second clustering data, the first neighborhood data, and the second neighborhood data.

2. The method according to claim 1, characterized in that: The first feature set includes first gene expression data and first site coordinate data, and the second feature set includes second gene expression data and second site coordinate data; The determining feature difference data according to the first feature set and the second feature set includes: determining a first gene expression difference according to the first gene expression data and the second gene expression data; Determine a first point distance difference according to the first point coordinate data and the second point coordinate data; Determine a first site weight matrix of the source slice and a second site weight matrix of the target slice respectively; the first site weight matrix is ​​used to describe the relative importance of each source site in the source slice; the second site weight matrix is ​​used to describe the relative importance of each target site in the target slice; Determining a first site weight difference according to the first site weight matrix and the second site weight matrix; The feature difference data is determined according to the first gene expression difference, the first site distance difference and the first site weight difference.

3. The method according to claim 1, characterized in that The first feature set includes first gene expression data and first site coordinate data, and the second feature set includes second gene expression data and second site coordinate data; The determining feature difference data according to the first feature set and the second feature set includes: Sampling the source site to obtain a source sampling site, and sampling the target site to obtain a target sampling site; Determine a second gene expression difference according to the source sampling site, the target sampling site, the first gene expression data, and the second gene expression data; Determine a second site distance difference according to the source sampling site, the target sampling site, the first site coordinate data, and the second site coordinate data; Determine a third site weight matrix of the source slice and a fourth site weight matrix of the target slice; the third site weight matrix is ​​used to describe the relative importance of each source site in the source slice; the fourth site weight matrix is ​​used to describe the relative importance of each target site in the target slice; Determining a second site weight difference according to the third site weight matrix and the fourth site weight matrix; The feature difference data is determined according to the second gene expression difference, the second site distance difference and the second site weight difference.

4. The method according to claim 2 or 3, characterized in that: Before determining feature difference data according to the first feature set and the second feature set, the method further includes performing dimensionality reduction processing on the first gene expression data, including: performing normalization processing on the first gene expression data and the second gene expression data respectively; Merging the normalized first gene expression data and the second gene expression data according to a common gene dimension to obtain gene merged data; Performing dimensionality reduction processing on the gene merged data to obtain total gene dimensionality reduction data; Sub-gene dimensionality reduction data of the first gene expression data are extracted from the total gene dimensionality reduction data based on the cell labels of the source slice, and the sub-gene dimensionality reduction data are used as the dimensionality reduction result of the first gene expression data.

5. The method according to claim 3, characterized in that: The first site coordinate data is used to describe the first initial coordinates of each source site of the source slice in the initial space, and the second site coordinate data is used to describe the second initial coordinates of each target site of the target slice in the initial space; Before determining feature difference data according to the first feature set and the second feature set, the method further includes performing dimensionality reduction processing on the first site coordinate data, including: Determine a site distance matrix according to the first initial coordinates and the second initial coordinates; Performing dimensionality reduction processing on the site distance matrix to obtain a first initial low-dimensional coordinate corresponding to each source site of the source slice in the low-dimensional space; Determine the nearest neighbor distance between the first initial coordinate of each of the source sites and the first initial low-dimensional coordinate, and determine a reference distance based on a plurality of the nearest neighbor distances; The low-dimensional space is scaled according to the reference distance, the first target low-dimensional coordinates of each source site are obtained according to the scaled low-dimensional space, and the dimensionality reduction result of the first site coordinate data is obtained according to the first target low-dimensional coordinates.

6. The method according to claim 5, characterized in that The sampling of the source site to obtain the source sampling site comprises: Determining a segmentation size according to the reference distance and a preset scaling factor; Segmenting the low-dimensional space according to the segmentation size to obtain segmentation units; The source sampling point is obtained by sampling the source point according to the segmentation unit.

7. The method according to claim 1, characterized in that The registering the source slice and the target slice according to the target transmission plan includes: Acquire first centroid data of the source site and second centroid data of the target site according to the target transmission plan; The registration coordinates are calculated according to the first centroid data and the second centroid data, and the source slice and the target slice are registered according to the registration coordinates.

8. A multi-slice registration device, characterized in that: The device comprises: a feature acquisition unit, configured to acquire a first feature set of the source slice based on the spatial transcription data of the source slice, and acquire a second feature set of the target slice based on the spatial transcription data of the target slice; a feature difference determining unit, configured to determine feature difference data according to the first feature set and the second feature set; a distribution data determination unit, configured to determine first distribution data and second distribution data, wherein the first distribution data is used to describe the total quality transmitted from the source site in the source slice to the target site in the target slice, and the second distribution data is used to describe the total quality received by the target site; a registration unit, configured to construct a target transmission plan according to the feature difference data and the distribution difference between the first distribution data and the second distribution data, and to register the source slice and the target slice according to the target transmission plan; The step of constructing a target transmission plan according to the characteristic difference data and the distribution difference between the first distribution data and the second distribution data includes: constructing an initial transmission plan according to the characteristic difference data and the distribution difference between the first distribution data and the second distribution data; Constructing site registration mask data according to the initial transmission plan, and constructing the target transmission plan according to the site registration mask data; wherein the site registration mask data is used to describe the sites that need to be considered in the fine mapping; The initial transmission plan is used to describe the quality of transmission from the source site to the target site; The constructing site registration mask data according to the initial transmission plan includes: Constructing a cluster set according to the maximum transmission quality in the initial transmission plan, constructing first cluster data of the source slice according to the relationship between the source sampling site and the cluster set, and constructing second cluster data of the target slice according to the relationship between the target sampling site and the cluster set; Constructing first neighborhood data of the source slice and second neighborhood data of the target slice according to the cluster set and a preset distance; The site registration mask data is constructed according to the first clustering data, the second clustering data, the first neighborhood data, and the second neighborhood data.

9. 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 multi-slice registration method according to any one of claims 1 to 7 is implemented.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the multi-slice registration method according to any one of claims 1 to 7 is implemented.

11. 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 multi-slice registration method according to any one of claims 1 to 7.

Images