A cell niche generation modeling and virtual perturbation method

Abstract

The application relates to the field of artificial intelligence algorithms and discloses a cell ecological niche generation modeling and virtual disturbance method, which comprises the following steps: step S1: data input and preprocessing, receiving single-cell or spatial omics data, performing quality control, normalization and logarithmic conversion processing on an expression matrix, and removing low-quality cells or features; step S2: constructing a cell neighborhood and extracting multi-view features, determining the neighborhood range of each cell, in the application, multi-view input is constructed from three dimensions of cell self molecular features, neighborhood molecular features and neighborhood cell composition features, through the collaborative architecture of a shared encoder and multiple view-specific decoders, the intrinsic heterogeneity of cells and the changes induced by the microenvironment are simultaneously captured in a unified latent representation space, the overall characterization of the functional heterogeneity of cell ecological niches and spatial organization characteristics is realized, and therefore more fine and biologically interpretable ecological niche analysis results are obtained.

Description

Technical Field

[0001] This invention relates to the field of artificial intelligence algorithms, and in particular to a method for cell niche generation modeling and virtual perturbation. Background Technology

[0002] Cellular niches are fundamental functional units that characterize the composition, spatial structure, and functional state of cells within a specific tissue region. They regulate key cellular behaviors such as growth, proliferation, differentiation, and apoptosis, and play a central role in various biological processes, including maintaining tissue homeostasis and shaping disease progression. With the rapid development of single-cell transcriptome sequencing and spatial transcriptome sequencing technologies, researchers are now able to simultaneously acquire molecular expression profiles and spatial location information of cells at single-cell resolution, providing a data foundation for accurately identifying diverse cellular niches in tissues and analyzing their dynamic changes.

[0003] However, existing computational methods still have significant limitations in cellular niche analysis. In static identification, methods such as NicheCompass, scNiche, STAGATE, and GraphST, while capable of niche clustering by integrating molecular expression and spatial information, are mostly limited to static identification under single conditions, lacking robust cross-sample integration capabilities, susceptible to batch effects, and lacking the ability to generate and model niche characteristics. In generative modeling and perturbation prediction, deep generative models such as scVI, scANVI, scGen, and Biolord, as well as pre-trained basic models such as Geneformer, scGPT, and Nicheformer, while showing potential in single-cell data integration and perturbation response prediction, primarily focus on individual cells, neglecting the regulatory influence of the cell's microenvironment, and their ability to perform cross-condition dynamic modeling and prediction at the cellular niche level remains limited. In the joint analysis of single-cell and spatial omics, existing methods such as NiCo and ENVI mainly focus on cell type annotation transfer or gene expression completion, lacking the ability to systematically reconstruct cellular niche information. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method for cell niche generation modeling and virtual perturbation, which solves the above problems.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for cell niche generation modeling and virtual perturbation, comprising the following steps: Step S1: Data input and preprocessing. Receive single-cell or spatial omics data, perform quality control, normalization, and logarithmic transformation on the expression matrix, and remove low-quality cells or features. Step S2: Construct cell neighborhoods and extract multi-view features, determine the neighborhood range of each cell, and extract multi-view features of each cell from the neighborhood range. The multi-view features include: cell-specific molecular features. Neighborhood molecular characteristics and the compositional characteristics of neighboring cells ; Step S3: Construct conditional variable embeddings. Preprocess cell-related biological covariates and project the preprocessed covariates into a continuous embedding space using a conditional encoder to obtain the conditional variable embeddings. ; Step S4: Establish a multi-view conditional variational autoencoder model, based on the cell's own molecular characteristics. and condition variable embedding As input, the latent representation is inferred through the encoder. The encoder reconstructs the features of each view using multiple view-specific decoders; the encoder employs variational distribution. Approximate posterior distribution and employ reparameterization techniques for sampling; each view decoder uses a latent representation. Embedded with condition variables The concatenation of the images is used as input to reconstruct the features of the corresponding views. Step S5: Model training and optimization. The Adam optimizer is used to train the model, and the optimization objective is the total loss function.

[0006] in, Represents the total loss function. This indicates a loss of cell expression reconstitution. The neighborhood represents the reconstruction loss. Indicates the reconstruction loss of the neighborhood composition. This represents the KL divergence loss. and Indicates the loss weight parameter; Step S6: Cellular niche identification and dynamic prediction, based on the latent representation obtained during training. Unsupervised clustering was performed to identify cell niches; and conditional variable embeddings were modified. Simultaneously maintain potential representations Without changing the input, the corresponding view decoder is used to generate counterfactual data, enabling dynamic prediction of cell niches, including prediction of cell niche feature expression, prediction of changes in cell niche cell neighborhood composition, and reconstruction of missing niche features.

[0007] According to one embodiment of the present invention, in step S2, the determination of the neighborhood range of each cell adopts the K nearest neighbor algorithm, and the number of nearest neighbor cells k is 20~40, with the optimum being 30; or a radius-based neighborhood definition strategy is adopted to include all cells within a preset spatial distance range into the neighborhood.

[0008] According to one embodiment of the present invention, in step S3, the biological covariate includes at least one of the following: batch, sample, age, and spatial background; for discrete covariates, one-hot encoding is used for binarization; for continuous covariates, z-score normalization or min-max scaling is used for normalization.

[0009] According to one embodiment of the present invention, in step S3, the conditional encoder is implemented using a linear layer or a multilayer perceptron, and the parameters ψ of the conditional encoder are jointly optimized with the parameters of the multi-view conditional variational autoencoder model.

[0010] According to one embodiment of the present invention, in step S4, the encoder includes three hidden layers with 2048, 1024 and 512 neurons respectively. Each hidden layer is followed by a ReLU activation function and a dropout layer with a dropout rate of 0.05~0.2, with an optimal rate of 0.1. The hidden layer structure of each view decoder is symmetrical to that of the encoder, with hidden layer dimensions of 512, 1024 and 2048 respectively.

[0011] According to one embodiment of the present invention, in step S4, the latent variable Obeying standard isotropic Gaussian priors The variational distribution approximates the posterior distribution as a multivariate Gaussian distribution with a diagonal covariance structure. Represents the identity matrix; the reparameterization technique is expressed as: ,in, Denotes the variational posterior distribution. Let these represent the mean and standard deviation of the latent variable, respectively. and By parameters The encoder network output, It represents the Hadamardi (or Hadama) stack.

[0012] According to one embodiment of the present invention, the neighborhood composition reconstruction loss adopts the cross-entropy loss function, and the neighborhood composition features are processed by the softmax function and then input into the cross-entropy loss calculation; the value range of the loss weight parameter α is 0.5~2, with the optimum being 1; the value range of the loss weight parameter β is 0.05~0.5, with the optimum being 0.1.

[0013] According to one embodiment of the present invention, in step S5, the number of epochs for model training is 500~1500, with an optimal value of 1000; the learning rate is 0.00005~0.0005, with an optimal value of 0.0001; and the batch size is 256~1024, with an optimal value of 512.

[0014] According to one embodiment of the present invention, in step S6, the unsupervised clustering adopts the k-means algorithm; when the target number of clusters is not specified, the optimal number of clusters is determined based on cluster stability, specifically: within the specified number of clusters, clustering is performed for each number of clusters K, the average Fowlkes-Mallows index between clusters K-1 and K and between clusters K and K+1 is calculated, and the number of clusters with the highest average Fowlkes-Mallows index is selected as the optimal number of clusters.

[0015] According to one embodiment of the present invention, the method supports joint modeling of single-cell omics data and spatial omics data. In joint modeling, a unified input space is constructed with shared genes, the encoder uses shared genes to learn cross-modal potential representations, the decoder reconstructs the feature information corresponding to spatial omics and single-cell omics respectively, and the technology platform information is used as a conditional variable input to the conditional encoder to eliminate systematic bias between different sequencing platforms.

[0016] Beneficial effects This invention provides a method for modeling and virtually perturbing cell niche generation. Compared with existing technologies, it has the following advantages: 1. In this invention, a conditional encoder is used to uniformly model multiple biological covariates such as batch, sample, age, and spatial background. This enables the model to learn latent representations that are insensitive to confounding covariates. While effectively eliminating batch effects, it accurately identifies conserved and condition-specific niches across samples, significantly improving the comparability and consistency of niche representations across different datasets. Furthermore, this method constructs a multi-view input from three dimensions: cell-specific molecular features, neighborhood molecular features, and neighborhood cell composition features. Through a collaborative architecture of a shared encoder and multiple view-specific decoders, it simultaneously captures changes in intrinsic cellular heterogeneity and microenvironment-induced variations within a unified latent representation space. This achieves a comprehensive characterization of cellular niche functional heterogeneity and spatial organization features, resulting in more refined and biologically interpretable niche analysis results. 2. This invention combines a niche-aware modeling strategy with a conditional variational autoencoder (CDAE) framework, overcoming the limitations of existing cell-centric perturbation models. It achieves dynamic prediction and counterfactual inference at the niche level. By modifying the embedding of conditional variables while maintaining the underlying representation, this method can predict changes in the characteristic expression of cell niches of interest under specific conditions, remodeling of cell neighborhood composition, and reconstruction of missing niche features. This extends traditional niche analysis from static identification to dynamic prediction and characterization. Simultaneously, this method supports joint modeling of single-cell omics and spatial omics data. By introducing technology platform information as a conditional variable to eliminate systematic bias, it significantly improves the analytical resolution of cell niches in complex tissues while achieving missing gene expression inference and reconstruction of missing niche information. This provides a powerful computational analysis tool for developmental biology, aging research, tumor immune microenvironment analysis, and precision treatment strategy formulation. Attached Figure Description

[0017] Figure 1 This is a flowchart of a cell niche generation modeling and virtual perturbation method proposed in this invention; Figure 2 Image showing the results of developmental brain fine-grained cell niche identification using nicheMAP; Figure 3 This is a diagram showing the prediction results of the dynamic changes in cellular niches during brain aging using nicheMAP. Figure 4 This is a diagram showing the prediction results of the dynamic changes in the tumor microenvironment T cell niche in nicheMAP. Figure 5 Figure showing the results of the joint modeling application of single-cell and spatial omics in nicheMAP; Figure 6 The figure shows the performance comparison results of nicheMAP and existing techniques on multiple spatial transcriptome datasets. Figure 7 The figure shows the performance comparison results of nicheMAP and the standard VAE model on multiple spatial transcriptome datasets. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] Please see Figures 1-7 The present invention provides the following technical solutions, specifically including the following embodiments: Example 1: A method for modeling and virtually perturbing cell niche generation includes the following steps: Step S1: Data input and preprocessing. Receive single-cell or spatial omics data, perform quality control, normalization, and logarithmic transformation on the expression matrix, and remove low-quality cells or features. Step S2: Construct cell neighborhoods and extract multi-view features. Determine the neighborhood range of each cell and extract multi-view features for each cell from the neighborhood range. Multi-view features include: cell-specific molecular features. Neighborhood molecular characteristics and the compositional characteristics of neighboring cells The neighborhood range of each cell is determined by the K nearest neighbor algorithm, where the number of nearest neighbor cells k is 20 to 40, with 30 being optimal; or a radius-based neighborhood definition strategy is adopted, which includes all cells within a preset spatial distance range into the neighborhood. Step S3: Construct conditional variable embeddings. Preprocess cell-related biological covariates and project the preprocessed covariates into a continuous embedding space using a conditional encoder to obtain the conditional variable embeddings. Biological covariates include at least one of batch, sample, age, and spatial context; for discrete covariates, one-hot encoding is used for binarization; for continuous covariates, z-score normalization or min-max scaling is used for normalization. The conditional encoder is implemented using a linear layer or multilayer perceptron. The parameters of the conditional encoder are... Joint optimization of parameters with the multi-view conditional variational autoencoder model; Step S4: Establish a multi-view conditional variational autoencoder model, based on the cell's own molecular characteristics. and condition variable embedding As input, the latent representation is inferred through the encoder. The encoder reconstructs the features of each view using multiple view-specific decoders; the encoder employs variational distribution. Approximate posterior distribution and employ reparameterization techniques for sampling; each view decoder uses a latent representation. Embedded with condition variables The concatenation of the inputs is used to reconstruct the features of the corresponding views. The encoder contains three hidden layers with 2048, 1024, and 512 neurons, respectively. Each hidden layer is followed by a ReLU activation function and a dropout layer with a dropout rate of 0.05~0.2, with an optimal rate of 0.1. The hidden layer structure of each view decoder is symmetrical to that of the encoder, with hidden layer dimensions of 512, 1024, and 2048, respectively. Latent variables... Obeying standard isotropic Gaussian priors The variational distribution approximates the posterior distribution as a multivariate Gaussian distribution with a diagonal covariance structure. Represents the identity matrix; the reparameterization technique is used to represent it as: ,in, Denotes the variational posterior distribution. Let these represent the mean and standard deviation of the latent variable, respectively. and By parameters The encoder network output, It represents the Hadamardi (or Hadama) stack; Step S5: Model training and optimization. The Adam optimizer is used to train the model, and the optimization objective is the total loss function.

[0020] in, Represents the total loss function. This indicates a loss of cell expression reconstitution. The neighborhood represents the reconstruction loss. Indicates the reconstruction loss of the neighborhood composition. This represents the KL divergence loss. and The loss weight parameters represent the parameters used in the neighborhood reconstruction loss calculation, where the neighborhood composition features are processed by the softmax function and then input into the cross-entropy loss function. The value range is 0.5 to 2, with an optimal value of 1; loss weight parameter The value range is 0.05~0.5, with the optimum being 0.1; the number of epochs for model training is 500~1500, with the optimum being 1000; the learning rate is 0.00005~0.0005, with the optimum being 0.0001; the batch size is 256~1024, with the optimum being 512. Step S6: Cellular niche identification and dynamic prediction, based on the latent representation obtained during training. Unsupervised clustering was performed to identify cell niches; and conditional variable embeddings were modified. Simultaneously maintain potential representations Without changing the input, the corresponding view decoder is used to generate counterfactual data, realizing dynamic prediction of cell niches, including prediction of cell niche feature expression, prediction of changes in cell neighborhood composition of cell niches, and reconstruction of missing niche features. Unsupervised clustering adopts the k-means algorithm. When the target number of clusters is not specified, the optimal number of clusters is determined based on cluster stability. Specifically, within the specified number of clusters, clustering is performed for each number of clusters K, and the average Fowlkes-Mallows index between clusters K-1 and K and between clusters K and K+1 is calculated. The cluster with the highest average Fowlkes-Mallows index is selected as the optimal number of clusters.

[0021] The method supports joint modeling of single-cell omics data and spatial omics data. In joint modeling, a unified input space is constructed with shared genes. The encoder uses shared genes to learn cross-modal potential representations. The decoder reconstructs the feature information corresponding to spatial omics and single-cell omics respectively. The technology platform information is used as a conditional variable input to the conditional encoder to eliminate systematic bias between different sequencing platforms.

[0022] Example 2: Application of nicheMAP in identifying fine-grained cellular niches in the developing brain; To verify the application capability of nicheMAP in identifying fine-grained cellular niches across samples, this embodiment uses the mouse embryonic organogenesis seqFISH spatial transcriptome dataset for validation. This dataset contains three tissue slices from different embryos, among which brain tissue development has a well-defined spatial anatomical structure and rich regional heterogeneity, making it suitable as a validation scenario for fine-grained cellular niche identification.

[0023] Based on the cell annotation information provided in the original study, 11,039 cells derived from brain tissue were screened from all cells. The data were then preprocessed, including gene expression matrix normalization, spatial coordinate extraction, and cell neighborhood construction. For each cell, its local microenvironment was determined using the K-nearest neighbor method, with a neighborhood size set to 30 nearest neighbor cells. Based on this, cell-specific molecular features, average neighborhood molecular features, and neighborhood cell composition features were extracted to construct multi-view input data for cell niches. The above data was used to train a nicheMAP model to obtain the potential niche representation for each cell. Furthermore, the Leiden clustering algorithm in the Scanpy package was used for unsupervised clustering analysis of the potential representations, with a clustering resolution set to 0.5, to identify fine-grained cell niches in the developing brain. The results are as follows Figure 2As shown, the results indicate that the latent representations learned by nicheMAP can effectively eliminate batch effects between different tissue sections. Furthermore, the identified cellular niches are highly consistent with the actual spatial anatomical structures, exhibiting a continuous spatial distribution along the forebrain, midbrain, and hindbrain developmental axes. Further differential expression analysis of the identified niches revealed that classic forebrain development regulatory genes such as Lhx2 and Pax6 are specifically highly expressed in forebrain-related niches; the midbrain-related niches show significantly high expression of the En1 gene, which has been shown to play a key regulatory role in the formation of midbrain dopaminergic precursor cells; and the hindbrain-related niches exhibit specific expression characteristics of Sfrp1 and multiple Hox family genes, reflecting significant region specificity during hindbrain segmental development. These results are consistent with existing biological research. To further evaluate the ability of nicheMAP to resolve fine-grained cellular niches, the crucial midbrain-hindbrain boundary (MHB) region during brain development was analyzed. This region, formed by the boundary between Otx2-highly expressed and Gbx2-highly expressed regions, is an important signal center and tissue boundary during embryonic brain development. Analysis results show that the latent representations learned by nicheMAP can effectively separate cell populations located in the midbrain and hindbrain, thereby accurately identifying the MHB niche. In contrast, existing techniques struggle to resolve this fine spatial boundary. These results demonstrate that nicheMAP can achieve high-resolution, fine-grained cellular niche identification and resolution.

[0024] Example 3: Application of nicheMAP in predicting dynamic changes in cellular niches during brain aging; To verify the ability of nicheMAP to model and predict the dynamic changes of cellular niches over time, this embodiment uses the MERFISH spatial transcriptome dataset of mouse brain aging for analysis. This dataset covers 20 age stages from 3.4 months to 34.5 months, containing 1,453,144 cells, and can systematically reflect the dynamic changes of different cell types and different brain region niches during brain aging.

[0025] First, the data underwent standard preprocessing, including expression matrix normalization, cell type annotation, and spatial coordinate information extraction. Then, cell neighborhoods were constructed using the nicheMAP workflow, and cell-specific molecular features, neighborhood molecular features, and neighborhood cell composition features were extracted. Age was used as a conditional variable input to the conditional encoder to establish a mapping relationship between age and cell niche status. During model training, data from 16 age groups were selected as the training set: 3.4, 4.3, 5.4, 6.6, 9.8, 12.9, 15.8, 18.8, 19.8, 21.4, 24.6, 26.7, 28.5, 30.9, 32.6, and 34.5 months. The remaining four age groups not used in training (3.8, 15.5, 23.5, and 33.2 months) served as the test set to verify the model's predictive ability under unknown conditions. After model training, 3.4-month-old brain tissue cells were used as the initial input to the model. By changing the age information in the condition variables, counterfactual predictions were made for the cellular niche characteristics of unseen age stages. The model simultaneously predicted cell expression profiles, neighborhood molecular features, and neighborhood composition features, thereby reconstructing the overall state of cellular niches at different age stages.

[0026] The results are as follows Figure 3As shown, the results indicate that the latent representations learned by nicheMAP not only maintain the biological differences among major cell types such as neurons, astrocytes, microglia, and oligodendrocytes, but also further preserve the spatial organization structure at different niche scales. Compared to directly performing dimensionality reduction analysis on the original expression profile, further comparison of the predicted results with actual test data revealed that nicheMAP can accurately predict transcriptomic changes in multiple unseen age stages. At the cell type, sub-brain region, and brain region levels, the predicted expression profiles and actual expression profiles showed high consistency. To further evaluate the spatial heterogeneity of aging processes in different brain region niches, a common aging score (CAS) for each cell was calculated using a publicly reported common aging gene feature set, and the aging rate of different brain region niches was further estimated. The results show that nicheMAP accurately reproduces the significant asynchrony between different brain region niches during brain aging. Specifically, the corpus callosum and anterior commissure region showed the fastest aging rate, the lateral ventricle region also showed a relatively rapid aging trend, while the striatum and cerebral cortex showed a relatively slow aging process. This result is highly consistent with spatial aging patterns observed in real-world data. Further analysis revealed that even cells of the same cell type exhibited significantly different aging characteristics in different brain niches. For example, astrocytes, oligodendrocytes, and microglia located in the corpus callosum region showed more significant aging-related transcriptomic alterations than in other brain regions. These results demonstrate that nicheMAP accurately captures the differentiated evolutionary trajectories of different cellular niches during aging, enabling accurate modeling and prediction of dynamic changes in cellular niches over time.

[0027] Example 4: Application of nicheMAP in predicting T cell niche spatial heterogeneity in the tumor microenvironment; To verify the ability of nicheMAP to model and predict the spatial dynamics of cell niches, this embodiment uses the human non-small cell lung cancer (NSCLC) CosMx SMI spatial transcriptome dataset for analysis. This dataset contains 8 tissue sections from 5 patients, with a total of 68,301 T cells, showing significant spatial heterogeneity in T cell niches.

[0028] First, all T cells were extracted based on the cell annotation results provided in the original study, and the data underwent standard preprocessing, including expression matrix normalization, cell type annotation organization, and spatial coordinate information extraction. Then, cell neighborhoods were constructed according to the nicheMAP workflow, and cell-specific molecular features, neighborhood molecular features, and neighborhood cell composition features were extracted. Furthermore, the average distance between each T cell and its 30 nearest tumor cells was calculated using spatial coordinate information, and this distance was input into the nicheMAP as a spatial conditional variable to describe the spatial relationship between the T cell's microenvironment and the tumor region. After model training, the T cell farthest from the tumor region was selected as the initial state. By gradually changing the spatial conditional variable, the migration process of T cells from the distal tumor region to the adjacent tumor region was simulated, thus realizing a virtual spatial perturbation experiment.

[0029] The results are as follows Figure 4As shown, the prediction results reveal that as T cells gradually approach the tumor region, their transcriptomic characteristics exhibit a continuous and coordinated trend of change. This process is highly consistent with the actual changes in T cells at different spatial locations in the original data, indicating that the spatial heterogeneity of T cell niches learned by nicheMAP has clear biological significance, rather than being a simple statistical fitting result. Simultaneously, nicheMAP can also synchronously predict the dynamic changes in the cellular neighborhood composition of T cell niches. The results show that as T cells gradually approach the tumor region, the proportion of immune cells in their neighborhood continuously decreases, while the proportion of tumor cells and stromal cells such as fibroblasts gradually increases. This pattern of change is highly consistent with the known tissue structure characteristics of the tumor invasion margin region. Further comparison of T cell states in tumor-adjacent niches and tumor-distal niches reveals that the two types of niches have significantly different transcriptomic characteristics. T cells in tumor-adjacent niches highly express matrix remodeling-related genes such as MMP1, MMP3, and MMP14, as well as chemokines such as CXCL9 and CXCL10. They also exhibit upregulation of immune regulation and exhaustion-related genes such as MX1, OAS1, IFI27, LAG3, LGALS3, and FOXP3, suggesting a gradual transition to an exhaustion state. T cells in tumor-distal niches primarily express naive and memory T cell marker genes such as IL7R, CCR7, and SELL, and are enriched in pathways related to T cell activation and homeostasis. To further validate the generalization ability of nicheMAP, the trained model was directly applied to another independent human high-grade serous carcinoma (HGSC) spatial transcriptome dataset. This dataset contains 97 tissue sections and 17,153 T cells. Using the same spatial perturbation strategy as before, nicheMAP accurately reproduced the spatial variations in T cell niche transcriptional expression and neighborhood composition. Statistical analysis showed that the nicheMAP prediction results achieved over 70% consistency between the two experiments, indicating that nicheMAP can reveal the spatial heterogeneity of conserved T cell niches across different cancers, rather than being applicable only to specific datasets. These results demonstrate that nicheMAP accurately captures the heterogeneity of T cell niches in different spatial locations within the tumor microenvironment, and can accurately model and predict the dynamic changes in cell niches in a spatial dimension.

[0030] Example 5: Application of nicheMAP in joint modeling of single-cell and spatial omics; To verify the application capability of nicheMAP in the joint analysis of single-cell and spatial omics, this embodiment uses 10X Xenium spatial transcriptome data from human breast cancer and corresponding matched single-cell RNA sequencing (scRNA-seq) data for analysis. The 10X Xenium data contains 158,030 cells and 303 genes; the matched scRNA-seq data contains 26,031 cells and 16,076 genes. Because the two technology platforms have different data characteristics—the former provides precise spatial information but has a limited number of genes detected, while the latter provides complete transcriptome information but lacks spatial location—they are suitable as a joint modeling validation scenario.

[0031] First, shared gene features were extracted from both datasets, and a unified input space was constructed according to the joint modeling strategy proposed by nicheMAP. The model encoder learned cross-modal latent representations using shared genes, while the decoder reconstructed the feature information corresponding to spatial omics and single-cell omics, respectively. During training, technology platform information was introduced as a conditional variable to eliminate systematic bias between different sequencing platforms. After model training, the inference ability of nicheMAP for missing gene expression information was evaluated. A five-fold cross-validation scheme was used, randomly dividing 303 target genes into 5 groups. Four groups were used for training each time, and the expression of the remaining group was predicted. The results showed that the spatial expression patterns predicted by nicheMAP were highly consistent with the actual expression patterns. Several representative genes were able to accurately recover their cell type-specific expression locations. For example, CEACAM6 was mainly located in the ductal carcinoma in situ region, FASN was mainly expressed in invasive tumor cells, POSTN was mainly distributed in the stromal cell region, IL7R was mainly distributed in the T cell region, and KRT15 was mainly expressed in the KRT15-positive myoepithelial cell population. Further validation using 10X Visium data from matched tissue sections revealed that nicheMAP could still accurately predict the spatial expression patterns of genes not included in the 303 gene detection panel, and the prediction results had higher spatial resolution than the 10X Visium technology.

[0032] The results are as follows Figure 5As shown, in addition to gene expression completion, nicheMAP can also infer missing spatial niche information in single-cell sequencing data. After using the trained model to predict single-cell data, multiple tumor-associated niche structures were successfully reconstructed, including different types of tumor-enriched niches and typical tertiary lymphoid structures (TLS), consistent with the actual spatial tissue structure. Further analysis revealed that nicheMAP can resolve compositional differences in different tumor-associated niches. For example, endothelial cells are widely distributed around various tumor niches, while ACTA2-positive and KRT15-positive myoepithelial cells are preferentially distributed around ductal carcinoma in situ niches. Furthermore, compared to DCIS#1, DCIS#2 niche showed more invasion-related cells, suggesting stronger invasive potential. Based on the whole transcriptome information inferred from nicheMAP, 4,110 endothelial cells located in different tumor niches were further selected for subpopulation analysis. The results identified nine endothelial cell subpopulations. This included classic capillary endothelial cells, venous endothelial cells, arterial endothelial cells, and lymphoid endothelial cells, as well as five tumor-associated endothelial cell subsets with significant tumor-related characteristics. Functional analysis showed that these subsets are involved in different biological processes such as angiogenesis, cell proliferation, extracellular matrix remodeling, antigen presentation, and immune activation. Particularly in invasive tumor niches, two endothelial cell subsets with angiogenesis and proliferation characteristics were significantly enriched and exhibited unique spatial distribution patterns. In contrast, while analysis using only raw 10X Xenium data could identify traditional vascular endothelial cell types, it was difficult to further analyze the aforementioned tumor-associated endothelial cell subsets. These results demonstrate that nicheMAP can achieve joint modeling of single-cell and spatial omics data, simultaneously performing multiple tasks such as inferring missing gene expression, reconstructing missing niche information, and resolving heterogeneity in complex niches, thereby significantly improving the resolution of cellular niche analysis in complex tissues.

[0033] Comparative Example 1: A comparative experiment between nicheMAP and existing technical methods; To verify the technical effectiveness of nicheMAP, several publicly available spatial transcriptome datasets with real cell niche annotation information were used for comparative testing, including the mouse primary visual cortex BaristaSeq dataset, the mouse medial prefrontal cortex STARmap dataset, the mouse hypothalamus anterior visual area MERFISH dataset, and the mouse whole brain STARmap PLUS dataset. Given the lack of methods specifically for cell niche generation modeling and virtual perturbation analysis in the existing technology, this comparative example selected three representative existing technical solutions as control groups: (1) deep generative models for single-cell generation modeling and perturbation prediction, including scVI, scANVI, and scGen; (2) unsupervised learning methods for static cell niche identification in spatial transcriptome data, including STAGATE and GraphST; and (3) basic pre-trained models that can be used for cell representation learning, generative modeling, or static niche identification tasks, including Geneformer, scGPT, and Nicheformer. All the methods described above are trained using their publicly recommended parameters and tested on the same dataset, with the same preprocessing procedures and evaluation system to ensure the fairness and reproducibility of the comparison results. To comprehensively evaluate the performance of different methods, they are assessed from the following four aspects: (1) Feature reconstruction accuracy; Feature reconstruction accuracy is used to evaluate the ability of the model to reconstruct cell and microenvironment features. Pearson correlation coefficient (PCC), cosine similarity (CS) and root mean square error (RMSE) were used for evaluation, and statistical analysis was performed at the gene level and the ecological niche level, respectively.

[0034] (2) Niche separation ability; Niche separation ability is used to evaluate whether the potential representations learned by the model can effectively distinguish different cell niches. The evaluation is based on indicators such as percentage of abnormal spots (PAS), spatial chaos score (CHAOS), average niche profile coefficient (ASW), graph cLISI, and graph connectivity.

[0035] (3) Niche biological information retention capacity; Niche biological information retention capacity is used to evaluate the consistency between potential representations and actual niche labels. Adjusted RAND Index (ARI), Normalized Mutual Information (NMI), Mutual Information (MI), and Logistic Score are used for evaluation. All potential representations learned by various methods are subjected to a unified clustering process, using the k-means algorithm for unsupervised clustering.

[0036] (4) Batch effect correction capability across samples; Batch effect correction capability is used to evaluate the model's ability to eliminate technical biases from different tissue sections and sample sources. The batch ASW, iLISI, and principal component regression (PCR) indicators were used for evaluation.

[0037] The results are as follows Figure 6 As shown, experimental results demonstrate that nicheMAP achieves the highest overall performance score across all test datasets and outperforms existing techniques in cell niche representation learning, feature reconstruction, and batch effect correction. Further analysis reveals that while existing cell niche identification methods can achieve niche clustering, they generally lack the ability to generate and model cell niche features, making it difficult to simultaneously address niche identification and feature reconstruction tasks. Furthermore, existing single-cell generation models, due to the lack of explicit modeling of spatial microenvironment information, cannot accurately capture tissue structure features at the niche level. Although the basic pre-trained model possesses some transfer learning capability, it still fails to reach the performance level of nicheMAP in the cell niche representation learning task. These results indicate that the multi-view conditional generation modeling framework proposed in this method can simultaneously learn the intrinsic state of cells and their microenvironment features, thereby significantly improving the accuracy of cell niche identification and the ability to reconstruct niche features, demonstrating superior technical performance compared to existing techniques.

[0038] Comparative Example 2: Ablation verification experiment of key technical features of nicheMAP; To verify the contribution of the multi-view niche modeling strategy in the nicheMAP method to the improvement of model performance, the following experimental scheme was constructed while maintaining consistency in training data, model parameters, and evaluation metrics: The control group used a standard variational autoencoder (VAE) structure, utilizing only the cell's own gene expression information as input and reconstruction target, without introducing the neighborhood molecular feature view or the neighborhood cell composition view; the experimental group used the nicheMAP model, simultaneously introducing the cell's own molecular features, neighborhood molecular features, and neighborhood cell composition features, and used a multi-view decoder for joint modeling. Both groups were trained and tested on the same dataset, and the same evaluation metrics as in Comparative Example 1 were used for performance evaluation. The results are as follows Figure 7As shown, experimental results indicate that the latent representations learned by the standard VAE model primarily reflect cell type differences but fail to effectively capture the structural information of the cell's microenvironment. Therefore, its performance in both cell niche identification and niche feature reconstruction tasks is significantly lower than that of nicheMAP. In contrast, nicheMAP, by simultaneously reconstructing multiple feature views related to the microenvironment, enables the model to learn the relationships between cells and their surrounding neighborhoods in the latent space. This results in a more realistic niche representation that better reflects the true tissue structure, significantly outperforming the standard VAE model in terms of niche separation ability and preservation of niche biological information. These results demonstrate that the multi-view niche modeling strategy proposed by nicheMAP is a key technical feature for achieving high-precision cell niche resolution and plays a significant role in improving model performance.

[0039] Conclusion: In summary, the cell niche generation modeling and virtual perturbation method proposed in this invention breaks through the bottleneck of existing technologies that are cell-centered and limited to static identification. It extends cell niche analysis from static characterization under a single condition to dynamic prediction and counterfactual generation across samples and conditions. This method has demonstrated excellent performance and wide applicability in multiple scenarios, including developmental biology, aging research, tumor immune microenvironment analysis, drug response prediction, and single-cell-spatial omics joint analysis. It provides a powerful computational tool and analytical framework for understanding the regulatory mechanisms of the cell microenvironment, revealing the laws of disease occurrence and development, and guiding precision treatment strategies.

[0040] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for modeling and virtually perturbing cell niche generation, characterized in that: Includes the following steps: Step S1: Data input and preprocessing. Receive single-cell or spatial omics data, perform quality control, normalization, and logarithmic transformation on the expression matrix, and remove low-quality cells or features. Step S2: Construct cell neighborhoods and extract multi-view features, determine the neighborhood range of each cell, and extract multi-view features of each cell from the neighborhood range. The multi-view features include: cell-specific molecular features. Neighborhood molecular characteristics and the compositional characteristics of neighboring cells ; Step S3: Construct conditional variable embeddings. Preprocess cell-related biological covariates and project the preprocessed covariates into a continuous embedding space using a conditional encoder to obtain the conditional variable embeddings. ; Step S4: Establish a multi-view conditional variational autoencoder model, based on the cell's own molecular characteristics. and condition variable embedding As input, the latent representation is inferred through the encoder. The encoder reconstructs the features of each view using multiple view-specific decoders; the encoder uses variational distribution. Approximate posterior distribution and employ reparameterization techniques for sampling; each view decoder uses a latent representation. Embedded with condition variables The concatenation of the images is used as input to reconstruct the features of the corresponding views. Step S5: Model training and optimization. The Adam optimizer is used to train the model, and the optimization objective is the total loss function. in, Represents the total loss function. This indicates a loss of cell expression reconstitution. The neighborhood represents the reconstruction loss. Indicates the reconstruction loss of the neighborhood composition. This represents the KL divergence loss. and Indicates the loss weight parameter; Step S6: Cellular niche identification and dynamic prediction, based on the latent representation obtained during training. Unsupervised clustering was performed to identify cell niches; and conditional variable embeddings were modified. Simultaneously maintain potential representations Without changing the input, the corresponding view decoder is used to generate counterfactual data, enabling dynamic prediction of cell niches, including prediction of cell niche feature expression, prediction of changes in cell niche cell neighborhood composition, and reconstruction of missing niche features.

2. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: In step S2, the neighborhood range of each cell is determined using the K nearest neighbor algorithm, where the number of nearest neighbor cells k is 20 to 40, with an optimal value of 30; or a radius-based neighborhood definition strategy is adopted to include all cells within a preset spatial distance range into the neighborhood.

3. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: In step S3, the biological covariates include at least one of batch, sample, age, and spatial background; for discrete covariates, one-hot encoding is used for binarization; for continuous covariates, z-score normalization or min-max scaling is used for normalization.

4. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: In step S3, the conditional encoder is implemented using a linear layer or a multilayer perceptron, and the parameters of the conditional encoder are... The parameters of the multi-view conditional variational autoencoder model are jointly optimized.

5. The cell niche generation modeling and virtual perturbation method according to claim 1, characterized in that: In step S4, the encoder contains three hidden layers with 2048, 1024 and 512 neurons respectively. Each hidden layer is followed by a ReLU activation function and a dropout layer with a dropout rate of 0.05~0.2, with an optimal rate of 0.

1. The hidden layer structure of each view decoder is symmetrical to that of the encoder, with hidden layer dimensions of 512, 1024 and 2048 respectively.

6. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: In step S4, the latent variables Obeying standard isotropic Gaussian priors The variational distribution approximates the posterior distribution as a multivariate Gaussian distribution with a diagonal covariance structure. Represents the identity matrix; the reparameterization technique is expressed as: ,in, Describes the variational posterior distribution. Let these represent the mean and standard deviation of the latent variable, respectively. and By parameters The encoder network output, It represents the Hadamardi (or Hadama) stack.

7. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: The neighborhood composition reconstruction loss uses the cross-entropy loss function, and the neighborhood composition features are processed by the softmax function and then input into the cross-entropy loss calculation; the loss weight parameters... The value range is 0.5 to 2, with an optimal value of 1; the loss weight parameter The value range is 0.05 to 0.5, with the optimal value being 0.

1.

8. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: In step S5, the number of epochs for model training is 500~1500, with an optimal value of 1000; the learning rate is 0.00005~0.0005, with an optimal value of 0.0001; and the batch size is 256~1024, with an optimal value of 512.

9. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: In step S6, the unsupervised clustering adopts the k-means algorithm. When the target number of clusters is not specified, the optimal number of clusters is determined based on the clustering stability. Specifically, within the specified number of clusters, clustering is performed for each number of clusters K, the average Fowlkes-Mallows index between clusters K-1 and K and between clusters K and K+1 is calculated, and the number of clusters with the highest average Fowlkes-Mallows index is selected as the optimal number of clusters.

10. The method for cell niche generation modeling and virtual perturbation according to claim 1, characterized in that: The method supports joint modeling of single-cell omics data and spatial omics data. In joint modeling, a unified input space is constructed using shared genes. The encoder uses shared genes to learn cross-modal potential representations. The decoder reconstructs the feature information corresponding to spatial omics and single-cell omics respectively. The technical platform information is used as a conditional variable input to the conditional encoder to eliminate systematic bias between different sequencing platforms.

Images