A method for spatial domain identification and batch effect removal in multi-slice spatial transcriptomics datasets

Abstract

The present invention relates to the technical field of cluster analysis of multi-slice spatial transcriptomics data sets, and discloses a method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics data sets. The present invention uses multi-layer stacking of graph neural networks as the backbone module of the feature encoder. In order to make the model more robust, the mask mechanism is integrated into the backbone module. At the same time, the present invention will also calculate the target function for those masked genes, and the normalized cosine error enhances the stability of embedded representation learning; secondly, in order to improve the spatial domain recognition performance of multi-slices, the present invention uses the anti-learning idea to construct a slice discriminator, which can effectively deal with the problem that multi-slice ST data cannot be accurately clustered; finally, the present invention also proposes a GGN greedy global information neighbor pair method to well alleviate the batch effect problem between slices caused by different technologies, different times and different conditions. The method of the present invention can not only effectively identify the spatial domain of multi-slice data, but also solve the batch effect problem between multi-slices.

Description

Technical Field

[0001] The present invention relates to the technical field of cluster analysis of multi-slice spatial transcriptomics datasets, and in particular to a method for spatial domain recognition and batch effect removal of multi-slice spatial transcriptomics datasets. Background Art

[0002] In complex organisms, cells are organized into similar clusters in a spatial environment. The complex arrangement of tissues reflects the dynamic interactions and special functions between cells, among which the precise location clues of transcriptional expression play a key role in revealing biological functions and networks. The latest spatial transcriptomics (Saptial Transcriptomics, ST) technology can provide dual functions. For example, 10x Visium achieves a spatial resolution of 55μm, and each capture site contains approximately 1-10 cells. It can perform whole genome analysis at the cellular and even subcellular level, while capturing gene expression corresponding to specific location points (called spots). These spatial insights provide a solid foundation for understanding many biological processes that affect health and disease, and have great application prospects in advancing diagnosis and treatment.

[0003] Identifying shared and specific spatial domains (clusters with similar spatial expression patterns) is one of the important tasks for comprehensive analysis of ST datasets. Currently, most methods are limited to clustering analysis of ST data on a single two-dimensional tissue slice. For example, K-means and Louvain algorithms only use gene expression data as input, which may lead to incoherent clustering results within tissue slices. SpaGCN uses gene expression, histological images and spatial location information as input, and uses an unsupervised clustering algorithm to detect different spatial expression patterns. However, batch effects between multiple ST slice data generated under different conditions, technologies or development stages may mask actual biological signals, limit our ability to understand biological processes occurring between multiple tissue slices, and hinder data integration, which poses a great challenge to existing single ST slice spatial domain identification methods.

[0004] To this end, some integrated analysis tools for multiple slices have emerged, such as STAligner, which achieves spatial perception of ST datasets by combining STAGATE and the mutual nearest neighbor MNN method based on latent space neighborhoods. It integrates multiple ST datasets by generating batch-corrected point embeddings, but it may mistakenly construct MNN pairs belonging to different spatial domains. Splane uses the cell type composition and spatial coordinate information of points obtained by Spoint prediction ST technology to jointly analyze multiple ST slices and decipher the complex interactions of gene expression in tissue structure, but this is a multi-stage processing method that first deconvolutes cell types from single-cell test datasets, and then uses its components to use graph convolutional networks and adversarial learning algorithms to identify spatial domains across multiple slices. STitch3D first uses the iterative closest point (ICP) algorithm or the PASTE algorithm to align multiple slices through spatial point registration, and then uses the distance between slices and the two-dimensional slice coordinates to construct a global three-dimensional spatial neighborhood graph. Finally, the gene expression and spatial information of multiple slices are mapped to the latent space through a deep graph autoencoder, and a deconvolution device is introduced to alleviate the batch effect. However, it relies on the alignment process of cross-slice coordinates, which limits its applicability to ST slices obtained under different conditions and in cases where there is only partial overlap.

[0005] Therefore, the prior art still needs to be improved and developed. Summary of the invention

[0006] In order to overcome the shortcomings of the prior art, the purpose of the present invention is to provide a method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics data set, aiming to solve the problem that the prior art cannot effectively identify the spatial domain of multi-slice data and cannot effectively solve the batch effect problem between multiple slices.

[0007] The first aspect of the present invention provides a method for spatial domain recognition and batch effect removal of a multi-slice spatial transcriptomics data set, comprising: preprocessing and data enhancement processing of the gene expression matrices of the acquired multiple spatial transcription slices to obtain a slice enhanced gene expression matrix; constructing a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice; constructing a neural network model, the neural network model comprising a feature encoder, a feature decoder, a slice discriminator and an anchor point neighbor selection module; inputting the enhanced gene expression matrix and the block diagonal adjacency matrix A into the feature encoder of the neural network model, and outputting the feature latent representation H of the spot; using the feature latent representation H of the spot as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, setting the learning objective functions of the generation stage and the adversarial stage to train the neural network model, and obtaining a trained neural network model; using the enhanced gene expression matrix of the spatial transcriptome slice to be analyzed and the block diagonal adjacency matrix A as the input of the feature encoder in the trained neural network model, and obtaining a feature latent representation with spatial domain recognition accuracy and batch effect removal.

[0008] Optionally, in a first implementation of the first aspect of the present invention, the gene expression matrices of the acquired multiple spatial transcription slices are preprocessed and data enhanced to obtain a slice enhanced gene expression matrix, including the steps of: connecting the gene expression matrices of the multiple spatial transcription slices along the dimension of spot, using the function provided by the SCANPY library to perform regularization processing on the selected highly variable genes to obtain a preprocessed gene expression matrix, denoted as X∈R Ns×2000 , where N s Represents the sum of the number of spots of multiple spatial transcription slices; randomly and uniformly samples the preprocessed gene expression matrix with probability p to obtain the indicator matrix M, and randomly masks part of the information of the preprocessed gene expression matrix to obtain the slice enhanced gene expression matrix

[0009] Optionally, in a second implementation of the first aspect of the present invention, a block diagonal adjacency matrix A is constructed based on the spatial position of the spot on each spatial transcription slice, including the steps of: calculating the two-dimensional Euclidean distance between spots on the spatial transcription slice, using the K nearest neighbor algorithm to ensure that each spot has k nearest neighbors; having the adjacency matrix A on the nth spatial transcription slice n , if the jth spot is a neighbor of the ith spot, then A nij =1, otherwise A nij =0, based on the adjacency matrix A n Build the block diagonal adjacency matrix A.

[0010] Optionally, in a third implementation of the first aspect of the present invention, the enhanced gene expression matrix and the block diagonal adjacency matrix A are input into the feature encoder of the neural network model, and the feature latent representation H of the spot is output, comprising the steps of: the feature encoder comprises several layers of graph convolution layers, and a convolution operation is performed on the input enhanced gene expression matrix and the block diagonal adjacency matrix A through the several layers of graph convolution layers, and the convolution operation expression is: Among them, H (l) is the feature representation of all spots in the current layer, H (l+1) It is the feature representation of all spots after a convolution operation. I is the identity matrix, D~ is the degree matrix of A~, σ is the activation function, W (l) is a trainable parameter matrix of the current layer; the several layers of graph convolutional layers perform convolution operations on the input enhanced gene expression matrix and the block diagonal adjacency matrix A through the feature encoder to obtain the feature latent representation H of the spot.

[0011] Optionally, in a fourth implementation of the first aspect of the present invention, the feature latent representation H of the spot is used as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, and the learning objective functions of the generation stage and the adversarial stage are set to train the neural network model to obtain a trained neural network model, including the steps of: after using the feature latent representation H of the spot as the input of the feature decoder, outputting the reconstructed gene expression matrix Z; after using the feature latent representation H of the spot as the input of the slice discriminator, outputting the label Y of the slice batch; after using the feature latent representation H of the spot as the input of the anchor point neighbor selection module, outputting the positive and negative anchor point pairs of the spot on the slice; setting the learning objective function of the generation stage to: L total =L class , where L class is the loss function of the first stage of the slice discriminator, Among them, θ represents the cross entropy loss function, and Y~ is the real slice ID; the learning objective function of the adversarial stage is set as: L total =μ1L sce +μ2L disc +μ3L tri , where L sce The objective loss function when the feature decoder reconstructs the gene expression matrix Z, L disc is the loss function of the second stage of the slice discriminator, L triThe anchor point neighbor selection module uses the ternary point loss as the target loss function for eliminating the batch effect, and μ1, μ2, and μ3 are all weight hyperparameter factors; after training the neural network model based on the learning objective functions of the generation stage and the adversarial stage, a trained neural network model is obtained.

[0012] Optionally, in a fifth implementation of the first aspect of the present invention, after the feature latent representation H of the spot is used as the input of the feature decoder, the reconstructed gene expression matrix Z is output, comprising the steps of: after the feature latent representation H of the spot is used as the input of the feature decoder, the reconstructed gene expression matrix Z is output, comprising the steps of: Calculate and obtain the reconstructed gene expression matrix Z, where Z (0) =H; the feature decoder uses the scaled cosine error as the objective function during the training process, expressed as: Among them, |V| represents the number of nodes in the mask set, M i Represents the mask vector of the i-th spot, Z i represents the features reconstructed by the model for the i-th spot, X i represents the initial unmasked feature of the i-th spot, T represents the transpose of the matrix, and γ represents the scaling factor.

[0013] Optionally, in a sixth implementation of the first aspect of the present invention, after the feature latent representation H of the spot is used as the input of the slice discriminator, the label Y of the slice batch is output, including the steps of: after the feature latent representation H of the spot is input into the slice discriminator, the operation process performed is as follows: Y (l+1) =σ(W (l) Y (l) +b (l) ), where Y (l) represents the features of the lth layer, b (l) represents the bias term of the lth layer, W (l) Represents the weight of the lth layer, σ represents the activation function, and the calculation result is processed by the Softmax function to output the label Y of the slice batch; the slice discriminator is divided into two stages during the training process, wherein the first stage is to allow the slice discriminator to learn the distinguishing features between batches from the feature latent representation H, and the second stage is to make the slice discriminator unable to distinguish which slice the spot comes from.

[0014] The second aspect of the present invention provides a device for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics data set, including: a data processing module, used to pre-process and perform data enhancement processing on the gene expression matrices of multiple spatial transcription slices to obtain a slice enhanced gene expression matrix; a matrix construction module, used to construct a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice; a model construction module, used to construct a neural network model, the neural network model includes a feature encoder, a feature decoder, a slice discriminator and an anchor point neighbor selection module; a feature extraction module, used to combine the enhanced gene expression matrix and the block diagonal adjacency matrix A. The matrix A is input into the feature encoder of the neural network model, and the feature latent representation H of the spot is output; the model training module is used to use the feature latent representation H of the spot as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, and the learning objective functions of the generation stage and the adversarial stage are set to train the neural network model to obtain the trained neural network model; the batch effect removal module is used to use the enhanced gene expression matrix of the spatial transcriptome slice to be analyzed and the block diagonal adjacency matrix A as the input of the feature encoder in the trained neural network model to obtain the feature latent representation with spatial domain recognition accuracy and batch effect removal.

[0015] The third aspect of the present invention provides a device for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset, comprising: a memory and at least one processor, the memory storing computer-readable instructions, the memory and the at least one processor being interconnected by lines; the at least one processor calling the computer-readable instructions in the memory so that the device for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset performs each step of the method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset as described above.

[0016] The fourth aspect of the present invention provides a computer-readable storage medium, which stores computer-readable instructions. When the computer-readable storage medium is run on a computer, it enables the computer to execute the various steps of the spatial domain identification and batch effect removal method of the multi-slice spatial transcriptomics data set as described above.

[0017] Beneficial effects: The present invention provides a method for spatial domain recognition and batch effect removal of multi-slice spatial transcriptomics data sets, with multi-layer stacking of graph neural networks as the backbone module of the feature encoder. In order to make the model more robust, the mask mechanism is integrated into the backbone module. At the same time, the present invention will also calculate the objective function for those masked genes, and the normalized cosine error enhances the stability of embedded representation learning; secondly, in order to improve the spatial domain recognition performance of multi-slices, the present invention uses the anti-learning idea to construct a slice discriminator, which can effectively deal with the problem of inaccurate clustering in multi-slice ST data; finally, the present invention also proposes a GGN greedy global information neighbor pair method to well alleviate the batch effect problem between slices caused by different technologies, different times and different conditions. The method of the present invention can not only effectively identify the spatial domain of multi-slice data, but also solve the batch effect problem between multi-slices. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 A flow chart of a method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset provided in an embodiment of the present invention.

[0019] Figure 2 Schematic diagram of the principle of the method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics datasets of the present invention.

[0020] Figure 3 A is a schematic diagram of DLPFC data; B is four consecutive slices including three donors; C is the comparison of our method with six baseline methods including SPIRIAL, STitch3D, Splane, STAligner, SEDR and DeepST in terms of ARI, Accuracy and Consistency indicators; D is the accuracy on four slices on donor 3 and its spatial domain recognition results.

[0021] Figure 4 The figure is a comparison result of the method of the present invention and six other baselines on the Local Inverse Simpson Index (LISI).

[0022] Figure 5 Integrated UMAP visualization on donor 3 (serial sections, 151673-151676) for different approaches, where colors indicate spatial domain type (lower layer) and batch (upper layer).

[0023] Figure 6 A schematic diagram of the structure of a device for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset provided in an embodiment of the present invention.

[0024] Figure 7A schematic diagram of the structure of a device for spatial domain recognition and batch effect removal of a multi-slice spatial transcriptomics dataset provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0025] The embodiment of the present invention provides a method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics data set. The terms "first", "second", "third", "fourth", etc. (if any) in the specification and claims of the present invention 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 described here can be implemented in an order other than that illustrated or described here. In addition, the terms "including" or "having" and any variations thereof 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.

[0026] See also Figure 1 , Figure 1 A flowchart of a preferred embodiment of a method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics data set provided by the present invention is shown in the figure, which includes the steps of:

[0027] S10, performing preprocessing and data enhancement processing on the acquired gene expression matrices of the multiple spatial transcription slices to obtain a slice-enhanced gene expression matrix;

[0028] S20, building a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice;

[0029] S30, constructing a neural network model, wherein the neural network model includes a feature encoder, a feature decoder, a slice discriminator, and an anchor point neighbor selection module;

[0030] S40, inputting the enhanced gene expression matrix and the block diagonal adjacency matrix A into the feature encoder of the neural network model, and outputting the feature latent representation H of the spot;

[0031] S50, using the feature latent representation H of the spot as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, setting the learning objective functions of the generation phase and the confrontation phase to train the neural network model, and obtaining a trained neural network model;

[0032] S60, using the enhanced gene expression matrix of the spatial transcriptome slice to be analyzed and the block diagonal adjacency matrix A as inputs of the feature encoder in the trained neural network model to obtain a feature latent representation with spatial domain recognition accuracy and batch effect removal.

[0033] Specifically, if Figure 2 As shown, the present invention uses multi-layer stacking of graph neural networks as the backbone module of the feature encoder. In order to make the model more robust, the mask mechanism is integrated into the backbone module. At the same time, the present invention will also calculate the target function for those masked genes, and the normalized cosine error enhances the stability of embedded representation learning; secondly, in order to improve the spatial domain recognition performance of multiple slices, the present invention uses the anti-learning idea to construct a slice discriminator, which can effectively deal with the problem of inaccurate clustering in multi-slice ST data; finally, the present invention also proposes a GGN greedy global information neighbor pair method to well alleviate the batch effect problem between slices caused by different technologies, different times and different conditions. The method of the present invention can not only effectively identify the spatial domain of multi-slice data, but also solve the batch effect problem between multi-slices.

[0034] In some embodiments, the gene expression matrices of the acquired multiple spatial transcription slices are preprocessed and data enhanced to obtain a slice enhanced gene expression matrix, which specifically includes the following steps: connecting the gene expression matrices of the multiple spatial transcription slices along the spot dimension, using the function provided by the SCANPY library to perform regularization processing on the selected highly variable genes, and obtaining a preprocessed gene expression matrix, denoted as X∈R Ns×2000 , where N s Represents the sum of the number of spots of multiple spatial transcription slices; In general autoencoders, if no additional constraints are imposed, they may suffer from the notorious "identity transformation" problem, that is, simply copying the input features to the output without learning useful information, and randomly masking part of the information of the gene expression matrix will effectively prevent such a simple mapping from occurring. Therefore, the present invention randomly and uniformly samples the preprocessed gene expression matrix with probability p to obtain the indicator matrix M, thereby obtaining the slice enhanced gene expression matrix

[0035] In some embodiments, a block diagonal adjacency matrix A is constructed based on the spatial position of the spot on each spatial transcription slice, specifically comprising the following steps: on the spatial transcription slice, the two-dimensional Euclidean distance between the spots is calculated, and the K nearest neighbor algorithm is used to ensure that each spot has k nearest neighbors; on the nth spatial transcription slice, there is an adjacency matrix A n , if the jth spot is a neighbor of the ith spot, then A nij =1, otherwise A nij=0, based on the adjacency matrix A n Construct a block diagonal adjacency matrix A, which is expressed as:

[0036] In some embodiments, the enhanced gene expression matrix and the block diagonal adjacency matrix A are input into the feature encoder of the neural network model, and the feature latent representation H of the spot is output, which specifically includes the steps of: the feature encoder includes several layers of graph convolution layers, and the input enhanced gene expression matrix and the block diagonal adjacency matrix A are convolved through the several layers of graph convolution layers, and the convolution operation expression is: Among them, H (l) is the feature representation of all spots in the current layer, H (l+1) It is the feature representation of all spots after a convolution operation. I is the identity matrix, yes The degree matrix, σ is the activation function, W (l) is a trainable parameter matrix of the current layer; the several layers of graph convolution layers obtain the feature latent representation H of the spot after the input enhanced gene expression matrix and the block diagonal adjacency matrix A are convolved by the feature encoder. Specifically, the graph convolution layers show a strong ability in processing and aggregating information from neighboring nodes. They are good at capturing the relationship with neighboring nodes and obtaining more comprehensive representation information, thereby promoting the integration of gene expression and spatial location information. In this embodiment, the enhanced gene expression matrix The diagonal adjacency matrix A of multiple slices is used as input to learn the feature latent representation H of the spot. The feature latent representation H learned in this embodiment has important uses. First, in the training stage, it will be used as the feature decoder input for reconstructing the original gene expression; secondly, it is also the input of the slice discriminator; after training to a certain stage, it will also be used as the input of the anchor point neighbor selection module (i.e., the greedy global neighbor anchor point selection method, GGN) described in the present invention to select the triple of points; when the training stage is completed, we use the feature latent representation H for downstream task analysis such as spatial domain recognition and visualization.

[0037] In some embodiments, the feature latent representation H of the spot is used as the input of the feature decoder, the slice discriminator, and the anchor point neighbor selection module respectively, and the learning objective functions of the generation phase and the adversarial phase are set to train the neural network model to obtain a trained neural network model, including the steps of: after using the feature latent representation H of the spot as the input of the feature decoder, outputting the reconstructed gene expression matrix Z; after using the feature latent representation H of the spot as the input of the slice discriminator, outputting the label Y of the slice batch; after using the feature latent representation H of the spot as the input of the anchor point neighbor selection module, outputting the positive and negative anchor point pairs of the spot on the slice; setting the learning objective function of the generation phase to: L total =L class , where L class is the loss function of the first stage of the slice discriminator, Among them, θ represents the cross entropy loss function, and is the real slice ID; the learning objective function of the adversarial stage is set as: L total =μ1L sce +μ2L disc +μ3L tri , where L sce The objective loss function when the feature decoder reconstructs the gene expression matrix Z, L disc is the loss function of the second stage of the slice discriminator, L tri The anchor point neighbor selection module uses the ternary point loss as the target loss function for eliminating the batch effect, and μ1, μ2, and μ3 are all weight hyperparameter factors; after training the neural network model based on the learning objective functions of the generation stage and the adversarial stage, a trained neural network model is obtained.

[0038] Specifically, the feature decoder is also composed of a multi-layer GCN backbone module, with the feature latent representation H as the input of this component, the purpose is to obtain the reconstructed gene expression matrix Z. The process is as follows:

[0039] Among them is Z (0) =H; when the reconstructed feature representation Z is obtained, it will be used to calculate the reconstruction loss; in order to reconstruct the mask feature from the given partially observed input feature, this embodiment uses the scaled cosine error (SCE) as the objective function, and the normalized cosine error enhances the stability of the embedded representation learning. Under the predefined scale factor γ, the similarity between the reconstructed prediction and the original input is calculated only on the masked gene, and the mathematical formula is as follows:

[0040] where γ is fixed to 2 throughout the experiment to reduce the weight of contributions from simple samples during training, and |V| represents the number of nodes in the mask set, Mi Represents the mask vector of the i-th spot, Z i represents the features reconstructed by the model for the i-th spot, X i represents the initial unmasked feature of the i-th spot, T represents the transpose of the matrix, and γ represents the scaling factor. In this embodiment, M is the mask matrix constructed when doing data enhancement. Here, if M ij = 1, indicating that the j-th column element of the i-th spot should be masked (so the enhanced input written before is Xm = X(1-M), which means when M ij = 0, the element is to be retained). Therefore, M i In fact, it is the mask vector, which is the mask vector of the i-th spot. Secondly, the loss in this embodiment is to calculate whether the masked features are reconstructed. Therefore, the loss is not calculated for the unmasked areas, and only the loss is calculated for the areas where M=1. i , is the feature reconstructed by the model for the i-th spot; X i It is the original feature that was not covered at the beginning; combined with the above, X i Multiply M i It is the value of the masked position. Then, because we need to calculate the similarity between each spot, we need to transpose this product. This T represents the transpose of the matrix.

[0041] In order to reduce the batch effect between samples, this embodiment designs a slice discriminator, which is composed of a multi-layer perceptron and finally outputs the label of the slice batch. Where N is the number of slices contained in the ST dataset. Specifically, it uses the feature latent representation H of the spot as the input of the slice discriminator, and its operation process is as follows: Y (l+1) =σ(W (l) Y (l) +b (l) ), where Y (l) represents the features of the lth layer, b (l) represents the bias term of the lth layer, W (l) represents the weight of the lth layer, σ represents the activation function, and finally the Softmax function is used to process the final output.

[0042] In this embodiment, the slice discriminator training process is divided into two stages of learning. The first stage allows the discriminator to learn the distinguishing features between batches from the feature latent representation H; the second stage makes the slice discriminator unable to distinguish which slice the spot comes from. Different from the standard generative adversarial model idea, the input of both stages here is the feature latent representation H, and only the learning target has changed. We regard the first stage as the generation stage, the purpose of which is to discover the distinguishing features between slices from the latent representation so as to improve the discriminant ability of the discriminator. When the slice discriminator can accurately determine which sample the feature latent representation H belongs to, it means that there are obvious distinguishing features between samples, but the batch effect will increase at this time; we regard the second stage as the adversarial stage, the purpose of which is to fool the discriminator. At this time, the discriminator weight will not be updated, but only used to determine which slice the spot belongs to.

[0043] Therefore, in order to train a good discriminator that can discover the distinguishing features between slices and retain the original structure, we design the loss function of the first stage as: In the formula, θ represents the cross entropy loss function, and Y~ is the real slice ID. In order to deceive the well-trained discriminator and make the slice sample ID and the low-dimensional latent representation indistinguishable, we design the loss function of the second stage as:

[0044] After the pre-training phase, the output feature latent representation H has certain semantic information, and its dimension is much lower than the input dimension, but it still has batch effect. Therefore, this embodiment also needs to use the GGN method to select a certain number of triples N at a certain step length t. tri , which is used to reduce the batch effect problem encountered when the model handles multiple slices.

[0045] Specifically, the goal of this embodiment is to force the encoder to learn batch correction embedding of multiple ST data through triple learning. This embodiment defines a point triple consisting of a target anchor point, a positive anchor point, and a negative anchor point. It should be noted that in order to avoid potential batch effects and omissions, the triple is defined in the embedding space rather than in the original gene expression space. The positive anchor pairs formed are used to overcome batch effects, while the negative anchor pairs guide alignment in a discriminative manner. Different from the MNN method used by STAligner, this embodiment improves the three problems faced by MNN in multi-ST slices and proposes a new greedy neighbor positive anchor selection method (GGN) that considers the global situation. First, MNN is forced to find mutually adjacent spot pairs between each pair of slices. However, in the actual multi-slice ST data set, not any two slices have the same spatial domain, which will mistakenly construct MNN pairs between different spatial domains. In addition, MNN requires the calculation of mutual proximity between two slices based on latent space features, and this requirement is relatively strong, and it is easy to ignore non-MNN pairs from the same spatial domain. Finally, when the number of slices in the ST dataset increases, the number of times to calculate MNN will increase in a flat manner. Second, this will limit the application of this method to more complex ST datasets.

[0046] To this end, the greedy global nearest positive anchor selection method (GGN) proposed in this embodiment will effectively alleviate the above challenges. More specifically, for the target slice S n , there are v i Belongs to S n , extract the feature latent representation H learned by the feature encoder i , and the positive anchor point is selected from all other slices, with v j Not S n , whose characteristic latent representation is H j For a given target slice, the node v to be queried i , calculate its cosine similarity with nodes on all other slices:

[0047]

[0048] Given the similarity, calculate each target node v j The node set of K-nearest neighbors of ni In order to obtain a more global set of positive anchor points, global semantic information is obtained through clustering, which can find nodes that share global semantic information with anchor points. The latent representation is divided into k clusters using kmeans, denoted as G = {G1, ..., G k}, and c(h i )∈{1,...,k}, represents the node v of the target slice i Characteristics of hi Then we can obtain the node set C that belongs to the same spatial domain as the target node from the remaining slices. ni = {v j |v j ∈G c(zi)}, as a similarity point that shares semantics with the target node globally.

[0049] Therefore, this embodiment obtains the intersection P of the nearest neighbor point in the target node latent space and the semantically similar node from a global perspective. ni =B ni ∩C ni , and record the intersection as the global positive anchor point set of the target node on other slices. According to our experiments, we found that this will lead to another problem, which may cause pattern collapse. This can be reduced by randomly selecting some nodes from the anchor point set. Thanks to the random selection of some anchor point sets, the entire anchor point set is no longer overly dependent on the original cluster distribution. The selection of negative anchor point sets is based on fast acquisition as the standard, so we can randomly select |P from the point set ni | points as the negative anchor point set of the target node, denoted as N ni Finally, this example uses triple point loss as the target loss function to eliminate batch effects:

[0050] Where a, p and n represent the indexes of anchor points, positive points and negative points respectively, T is the set of triplets of target anchor points and positive and negative anchor points we selected, and N tri is the number of all triplets and τ (default 1.0) is the margin used to enforce distance between positive and negative pairs.

[0051] In this embodiment, the learning process of the neural network model is divided into two stages. The first is the pre-training stage. Since the model has not yet learned meaningful feature latent representations, calculating the neighbor pairs in the latent space at this time will make the model difficult to learn. Therefore, in the pre-training stage, the spatial domain recognition accuracy of the model is improved mainly through reconstruction loss and adversarial learning between multiple slices. Specifically, the following learning objectives are set in the generation stage of pre-training: L total =L class , at this time only the discriminator weight parameters are updated, and the pre-training adversarial stage has the following learning objectives: L total =μ1L sce +μ2L disc , where μ1, μ2 are the loss weight hyperparameter factors. At this time, the model updates the feature encoder and decoder weights without updating the discriminator weights.

[0052] After the pre-training phase, the latent space has certain semantic information and its dimension is much lower than the input dimension. We use the GGN method to select a certain number of triplets N at a certain step length t. tri , which is used to reduce the batch effect problem encountered when the model processes multiple slices. At this time, the learning objective of the generation phase does not change, while the learning objective of the adversarial phase is updated as follows: L total =μ1L sce +μ2L disc +μ3L tri , where μ3 is the weight hyperparameter factor of the triplet loss. When the training phase is over, the obtained feature latent space representation will be used for downstream task analysis.

[0053] In summary, the present invention uses multi-layer stacking of graph neural networks as the backbone module of the feature encoder. In order to make the model more robust, the mask mechanism is integrated into the backbone module. At the same time, the present invention will also calculate the target function for those masked genes, and the normalized cosine error enhances the stability of embedded representation learning; secondly, in order to improve the spatial domain recognition performance of multi-slices, the present invention uses the anti-learning idea to construct a slice discriminator, which can effectively deal with the problem of inaccurate clustering in multi-slice ST data; finally, the present invention also proposes a GGN greedy global information neighbor pair method to well alleviate the batch effect problem between slices caused by different technologies, different times and different conditions. The method of the present invention can not only effectively identify the spatial domain of multi-slice data, but also solve the batch effect problem between multi-slices.

[0054] The following specific experiments compare the capabilities and robustness of the method of the present invention (OURS) with other mainstream methods (SPIRIAL, STitch3D, Splane, STAligner, SEDR, DeepST). This example is analyzed on a 10x Vision dataset consisting of 12 slices of the human dorsolateral prefrontal cortex (DLPFC) containing spatial expression information. Maynard et al. manually annotated the layers of DLPFC and white matter (WM) based on morphological features and genetic markers, with 5-7 regions manually annotated for each slice.

[0055] Usually, indicators such as ARI (Adjusted Rand Index), NMI (Normalized Mutual Information), HOM (homogeneity), COM (completeness), CHAOS (spatial chaos score), PAS (percentage of outliers), and LISI (local inverse Simpson index) can be used to evaluate model performance.

[0056] 1. ARI (Adjusted Rand Index) is an indicator used to compare the similarity between two clustering results, taking into account the consistency between paired samples in the clustering results. Its calculation formula is as follows:

[0057] Among them, a i represents the number of sample pairs that belong to the same cluster in the i-th clustering result, b j represents the number of sample pairs that belong to the same category in the jth true label, n ij Indicates the number of sample pairs that are simultaneously assigned to the i-th clustering result and the j-th true label. It represents the number of combinations of selecting k elements from n elements, where N is the total number of samples. The value range of this indicator is between 0 and 1. The higher the value, the higher the similarity between the clustering result and the true label.

[0058] 2. NMI (Normalized Mutual Information) is an indicator used to evaluate clustering quality. It is based on the concept of information theory and measures the normalized value of the mutual information between the clustering result and the true label. Its calculation formula is as follows:

[0059] Among them, C represents the cluster set of clustering results, |C| represents the number of clusters, K represents the number of categories of true labels, and n ij represents the number of samples that belong to both cluster j and true label i, n i . represents the number of samples in cluster i, n. j represents the number of samples in the true label j, and N represents the total number of samples.

[0060] 3. The HOM (homogeneity) score is a metric that quantifies the homogeneity of clustered labeling when compared to known ground truth. A clustering result is considered homogeneous if all clusters contain only data points belonging to a single class. The HOM score is expressed as a value between 0 and 1, with 1 representing completely homogeneous labeling.

[0061] 4. The COM (completeness) score is a metric used to evaluate the completeness of the cluster labeling relative to the ground truth. If all data points belonging to a class are grouped into the same class, the clustering result is considered complete. The COM score ranges from 0 to 1, with 1 indicating perfect labeling.

[0062] We define here the accuracy index (Accuracy), which is the overall score based on the obtained NMI score, HOM score and COM score, and its calculation method is as follows:

[0063] The accuracy ranges from 0 to 1. When it approaches 1, it means that the unsupervised method can effectively cluster features, otherwise the effect is poor.

[0064] 5. The CHAOS (spatial chaos score) score is a metric used to evaluate the spatial continuity performance in the field of mass spectrometry imaging and spatial transcriptomics. The lower the chaos value, the better the continuity of spatial domain recognition. In order to apply CHAOS to quantify the performance of each spatial clustering method, we first construct a 1-nearest neighbor (NN) graph for each data. Specifically, each spot is connected to another spot with the minimum Euclidean distance in physical space. For this 1-NN graph, assuming d ij It is a cell in physical space i and cell j The Euclidean distance between , we calculate w as follows:

[0065] Assume n k is the number of cells in the kth spatial domain, n is the total number of cells in the data, and k is the number of unique spatial domains. CHAOS is calculated as:

[0066] 6. The PAS (Percentage of Outliers) score measures the randomness of points that are outside the spatial domain to which they are clustered. The PAS score is calculated as the proportion of points that have a different cluster label than at least 6 of the 10 adjacent points. The lower the PAS score, the better the continuity of the detected spatial domain, which means that the cells in the spatial domain are more homogeneous. Therefore, we have the calculation of the consistency composite score as follows, the lower the better the continuity:

[0067] 7. The LISI (Local Inverse Simpson Index) index can be used to evaluate the mixing of batches and spatial domains across multiple slices. LISI selects the nearest neighbors with a fixed perplexity based on the local distribution of distances. The selected neighbors are then used to calculate the inverse Simpson index of diversity, which is the effective number of types present in that neighborhood. The LISI (Local Inverse Simpson Index) of the batch in each domain is used to calculate the perplexity of the batch. batch ) is calculated for the batch labels, and a fraction close to the expected batch size indicates good mixing, as well as the Local Inverse Simpson Index (LISI) of the domain across all data domain ) A score close to 1 indicates that the cluster contains pure spatial domains to assess the degree of separation between the same domain and different domains. Normalized LISI batch and LISI domain The combination is calculated as: The higher the score, the better the batch correction effect.

[0068] Figure 3A is a schematic diagram of DLPFC data; B is four consecutive slices including three donors; C is the comparison of our method with six baseline methods including SPIRIAL, STitch3D, Splane, STAligner, SEDR and DeepST in terms of ARI, Accuracy and Consistency indicators; D is the accuracy on four slices on donor 3 and its spatial domain recognition results.

[0069] Since the DLPFC dataset consists of three tissue samples, each of which contains four slices, we take four slices as input at a time and finally count the accuracy and continuity of all 12 slices. Figure 3 The results show the accuracy and consistency of our spatial domain recognition on the DLPFC dataset. In the experimental results, we can see that our method achieves the best clustering accuracy and the best consistency. From the spatial domain recognition result diagram, we can see that each spatial domain has a very clear and continuous hierarchical structure. The SPIRIAL method has a mixture of spots in each domain, STitch3D incorrectly recognizes multiple domains, and STAligner cannot recognize Layer 5 and Layer 6.

[0070] We then compared the proposed method with six other baselines on the Local Inverse Simpson Index (LISI) (results are shown in Figure 4 ), and also visualized using UMAP integrated on donor 3 (consecutive sections, 151673-151676) using different methods (results shown in Figure 5 shown), in Figure 4 and Figure 5 The results of eliminating batch effects of the present invention are shown in the figure. It can be seen from the figure that our method has the best local inverse Simpson index median, and the batch mixing is better in UMAP, and the layer is very clearly identified. This experiment proves the effectiveness of our method. The above describes the spatial domain recognition and batch effect removal method of the multi-slice spatial transcriptomics dataset in the embodiment of the present invention. The following describes the spatial domain recognition and batch effect removal device of the multi-slice spatial transcriptomics dataset in the embodiment of the present invention. Please refer to Figure 6 In one embodiment of the present invention, a device for spatial domain recognition and batch effect removal of a multi-slice spatial transcriptomics data set includes:

[0071] The data processing module 10 is used to perform preprocessing and data enhancement processing on the gene expression matrices of the acquired multiple spatial transcription slices to obtain a slice enhanced gene expression matrix;

[0072] A matrix construction module 20, for constructing a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice;

[0073] A model building module 30, used to build a neural network model, wherein the neural network model includes a feature encoder, a feature decoder, a slice discriminator, and an anchor point neighbor selection module;

[0074] A feature extraction module 40, used for inputting the enhanced gene expression matrix and the block diagonal adjacency matrix A into a feature encoder of the neural network model, and outputting a feature latent representation H of a spot;

[0075] A model training module 50 is used to use the feature latent representation H of the spot as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module, respectively, and set the learning objective function of the generation phase and the adversarial phase to train the neural network model to obtain a trained neural network model;

[0076] The batch effect removal module 60 is used to use the enhanced gene expression matrix and the block diagonal adjacency matrix A of the spatial transcriptome slice to be analyzed as inputs of the feature encoder in the trained neural network model to obtain a feature latent representation with spatial domain recognition accuracy and batch effect removal.

[0077] above Figure 6 From the perspective of modular functional entities, the spatial domain identification and batch effect removal device of the multi-slice spatial transcriptomics dataset in the embodiment of the present invention is described in detail. The spatial domain identification and batch effect removal device of the multi-slice spatial transcriptomics dataset in the embodiment of the present invention is described in detail from the perspective of hardware processing.

[0078] Figure 7It is a structural schematic diagram of a spatial domain identification and batch effect removal device for a multi-slice spatial transcriptomics data set provided by an embodiment of the present invention. The spatial domain identification and batch effect removal device 100 for the multi-slice spatial transcriptomics data set may have relatively large differences due to different configurations or performances, and may include one or more processors (central processing units, CPU) 11 (for example, one or more processors) and a memory 12, and one or more storage media 13 (for example, one or more massive storage devices) storing application programs 133 or data 132. Among them, the memory 12 and the storage medium 13 can be short-term storage or permanent storage. The program stored in the storage medium 13 may include one or more modules (not shown in the figure), and each module may include a series of instruction operations in the spatial domain identification and batch effect removal device 100 for multi-slice spatial transcriptomics data sets. Furthermore, the processor 11 can be configured to communicate with the storage medium 13 to execute a series of instruction operations in the storage medium 13 on the prediction device 100 of the super-resolution gene expression map.

[0079] The spatial domain identification and batch effect removal device 100 for multi-slice spatial transcriptomics data sets may also include one or more power supplies 14, one or more wired or wireless network interfaces 15, one or more input and output interfaces 16, and / or one or more operating systems 131, such as Windows Serve, Mac OS X, Unix, Linux, FreeBSD, etc. It can be understood by those skilled in the art that Figure 7 The device structure shown does not constitute a limitation on the device 100 for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics datasets, and may include more or fewer components than shown in the figure, or combine certain components, or arrange the components differently.

[0080] The present invention also provides a computer-readable storage medium, which may be a non-volatile computer-readable storage medium or a volatile computer-readable storage medium. The computer-readable storage medium stores instructions, which, when executed on a computer, enable the computer to execute the steps of a method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics data set.

[0081] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-described system, device, or unit can refer to the corresponding process in the aforementioned method embodiment and will not be repeated here.

[0082] 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 invention is essentially or the part that contributes to the prior art or all or part of the technical solution can be embodied in the form of a software product. The computer software product is stored in a storage medium, including a number of instructions to enable a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM), random access memory (RAM), disk or optical disk and other media that can store program codes.

[0083] As described above, the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features thereof may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics datasets, characterized in that: Includes steps: Preprocessing and data enhancement processing are performed on the gene expression matrices of multiple spatial transcription slices to obtain a slice-enhanced gene expression matrix; Build a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice; Constructing a neural network model, wherein the neural network model includes a feature encoder, a feature decoder, a slice discriminator, and an anchor point neighbor selection module; Inputting the enhanced gene expression matrix and the block diagonal adjacency matrix A into the feature encoder of the neural network model, and outputting the feature latent representation H of the spot; The feature latent representation H of the spot is used as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, and the learning objective functions of the generation phase and the adversarial phase are set to train the neural network model to obtain a trained neural network model; The enhanced gene expression matrix and the block diagonal adjacency matrix A of the spatial transcriptome slice to be analyzed are used as inputs of the feature encoder in the trained neural network model to obtain a feature latent representation with spatial domain recognition accuracy and batch effect removal.

2. The method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics dataset according to claim 1, characterized in that: The gene expression matrices of the acquired multiple spatial transcription slices are preprocessed and data enhanced to obtain a slice enhanced gene expression matrix, including the following steps: The gene expression matrices of multiple spatial transcription slices are connected along the spot dimension, and the functions provided by the SCANPY library are used to regularize the selected highly variable genes to obtain the preprocessed gene expression matrix, denoted as X∈R Ns×2000 , where N s Represents the total number of spots in multiple spatial transcription slices; The pre-processed gene expression matrix is ​​randomly and uniformly sampled with probability p to obtain an indicator matrix M, and part of the information of the pre-processed gene expression matrix is ​​randomly masked to obtain a slice enhanced gene expression matrix 3. The method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics dataset according to claim 1, characterized in that: Building a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice includes the following steps: On the spatial transcription slice, the two-dimensional Euclidean distance between spots is calculated, and the K nearest neighbor algorithm is used to ensure that each spot has k nearest neighbors; On the nth spatial transcription slice there is an adjacency matrix A n , if the jth spot is a neighbor of the ith spot, then A nij =1, otherwise A nij =0, based on the adjacency matrix A n Build the block diagonal adjacency matrix A.

4. The method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics dataset according to claim 1, characterized in that: Inputting the enhanced gene expression matrix and the block diagonal adjacency matrix A into the feature encoder of the neural network model, and outputting the feature latent representation H of the spot, comprises the steps of: The feature encoder includes several graph convolution layers, through which the input enhanced gene expression matrix and the block diagonal adjacency matrix A are convolved. The convolution operation expression is: Among them, H (l) is the feature representation of all spots in the current layer, H (l+1) It is the feature representation of all spots after a convolution operation. I is the identity matrix, yes The degree matrix, σ is the activation function, W (l) is the trainable parameter matrix of the current layer; The plurality of graph convolution layers perform convolution operations on the input enhanced gene expression matrix and the block diagonal adjacency matrix A through a feature encoder to obtain a feature latent representation H of the spot.

5. The method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics dataset according to claim 1, characterized in that: The feature latent representation H of the spot is used as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, and the learning objective functions of the generation phase and the confrontation phase are set to train the neural network model to obtain a trained neural network model, including the steps of: After taking the feature latent representation H of the spot as the input of the feature decoder, the reconstructed gene expression matrix Z is output; After taking the feature latent representation H of the spot as the input of the slice discriminator, the label Y of the slice batch is output; After taking the feature latent representation H of the spot as the input of the anchor point neighbor selection module, the positive and negative anchor point pairs of the spot on the slice are output; The learning objective function of the generation phase is set as: L total =L class , where L class is the loss function of the first stage of the slice discriminator, Among them, θ represents the cross entropy loss function, and It is the real slice ID; The learning objective function of the adversarial stage is set as: L total =μ1L sce +μ2L disc +μ3L tri , where L sce The objective loss function when the feature decoder reconstructs the gene expression matrix Z, L disc is the loss function of the second stage of the slice discriminator, L tri The anchor point neighbor selection module uses the ternary point loss as the target loss function to eliminate the batch effect, and μ1, μ2, and μ3 are all weight hyperparameter factors; After training the neural network model based on the learning objective functions of the generation phase and the confrontation phase, a trained neural network model is obtained.

6. The method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics dataset according to claim 5, characterized in that: After taking the feature latent representation H of the spot as the input of the feature decoder, the reconstructed gene expression matrix Z is output, including the steps of: After taking the feature latent representation H of the spot as the input of the feature decoder, the formula Calculate and obtain the reconstructed gene expression matrix Z, where Z (0) =H; The feature decoder uses the scaled cosine error as the objective function during training, expressed as: Among them, |V| represents the number of nodes in the mask set, M i Represents the mask vector of the i-th spot, Z i represents the features reconstructed by the model for the i-th spot, X i represents the initial unmasked feature of the i-th spot, T represents the transpose of the matrix, and γ represents the scaling factor.

7. The method for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics dataset according to claim 5, characterized in that: After the feature latent representation H of the spot is used as the input of the slice discriminator, the label Y of the slice batch is output, including the steps: After the feature latent representation H of the spot is input into the slice discriminator, the operation process is as follows: (l+1) =σ(W (l) Y (l) +b (l) ), where Y (l) represents the features of the lth layer, b (l) represents the bias term of the lth layer, W (l) represents the weight of the lth layer, σ represents the activation function, the calculation result is processed by the Softmax function, and the label Y of the slice batch is output; The slice discriminator is divided into two stages during the training process, wherein the first stage is to allow the slice discriminator to learn the distinguishing features between batches from the feature latent representation H, and the second stage is to make the slice discriminator unable to distinguish which slice the spot comes from.

8. A device for spatial domain identification and batch effect removal of multi-slice spatial transcriptomics datasets, characterized in that: include: A data processing module is used to perform preprocessing and data enhancement processing on the gene expression matrices of the acquired multiple spatial transcription slices to obtain a slice-enhanced gene expression matrix; A matrix building module for building a block diagonal adjacency matrix A based on the spatial position of the spot on each spatial transcription slice; A model building module, used to build a neural network model, wherein the neural network model includes a feature encoder, a feature decoder, a slice discriminator, and an anchor point neighbor selection module; A feature extraction module, used for inputting the enhanced gene expression matrix and the block diagonal adjacency matrix A into the feature encoder of the neural network model, and outputting a feature latent representation H of a spot; A model training module, used to use the feature latent representation H of the spot as the input of the feature decoder, the slice discriminator and the anchor point neighbor selection module respectively, set the learning objective function of the generation phase and the adversarial phase to train the neural network model, and obtain a trained neural network model; The batch effect removal module is used to use the enhanced gene expression matrix and the block diagonal adjacency matrix A of the spatial transcriptome slice to be analyzed as the input of the feature encoder in the trained neural network model to obtain a feature latent representation with spatial domain recognition accuracy and batch effect removal.

9. A device for spatial domain recognition and batch effect removal of multi-slice spatial transcriptomics datasets, characterized in that: comprising a memory and at least one processor, wherein the memory has computer-readable instructions stored therein; The at least one processor calls the computer-readable instructions in the memory to execute each step of the method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset as described in any one of claims 1 to 7.

10. A computer-readable storage medium having computer-readable instructions stored thereon, characterized in that: When the computer-readable instructions are executed by a processor, the steps of the method for spatial domain identification and batch effect removal of a multi-slice spatial transcriptomics dataset are implemented.

Images