Method, system and computer device for screening of drug targets with spatiotemporal specificity

Abstract

The application discloses a screening method and system of a drug target with space-time specificity and a computer device, and relates to the technical field of bioinformatics and computational biology. The screening method is based on single-cell transcriptome sequencing data, and the method comprises the following steps: (1) quantitatively reconstructing the spatial positioning and functional mode of cells in a tissue, which comprises the following steps: 1.1) performing data preprocessing on the single-cell transcriptome sequencing data; 1.2) reconstructing the spatial positioning of single cells; and 1.3) reconstructing the biological functional mode of single cells; and (2) screening a drug target with space-time specificity, which comprises the following steps: 2.1) analyzing cell-cell communication; 2.2) constructing a gene regulatory network (GRN) with a specific tissue microenvironment state as the core; and 2.3) discovering a target. The application verifies an innovative single-cell data analysis method, and provides a new thought and technical path for drug research and development of metabolic diseases and other systemic diseases, and has the basis of popularization and application.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and computational biology, and in particular to a method, system, and computer device for screening drug targets with spatiotemporal specificity. Background Technology

[0002] Systemic diseases are diseases affecting multiple organs or systems throughout the body. Their causes often involve immune abnormalities, metabolic disorders, genetic factors, or infections, requiring comprehensive diagnosis and treatment. These diseases are usually not attributable to a single cause and have complex symptoms, necessitating multidisciplinary intervention. Systemic diseases include, for example, chronic metabolic diseases and tumors. Chronic metabolic diseases are a class of chronic illnesses caused by long-term abnormalities in metabolic processes, involving disorders of carbohydrate, fat, and protein metabolism, and are usually related to genetic factors and unhealthy lifestyles. Common types include diabetes, hypertension, hyperlipidemia, hyperuricemia (gout), and obesity. These diseases have an insidious onset, a long course, and can cause multi-organ damage, requiring long-term intervention and control.

[0003] One of the fundamental reasons why systemic diseases, such as chronic metabolic diseases and tumors, are difficult to cure is that they do not originate from a single cell type or molecular event, but are driven by the spatiotemporal interactions of multiple cell populations within a specific tissue microenvironment. The interwoven signaling networks of metabolism, immunity, and environmental stress within tissues, along with the spatiotemporal dynamics of tissue pathology, make the pathophysiological mechanisms of chronic metabolic diseases such as fatty liver disease extremely complex. For example, spatial location information is a crucial dimension for reconstructing the phenotype of the tissue microenvironment. The functional state of cells within a tissue is closely related to their spatial region, especially in organs with spatially partitioned structures. Cells in different spatial regions may be exposed to different blood flow components, oxygen partial pressures, inflammatory factor concentrations, and matrix microenvironments, thus exhibiting highly heterogeneous functional characteristics.

[0004] Therefore, developing spatiotemporally specific drug targets targeting specific tissue microenvironments, and addressing the stage-specific characteristics and spatial distribution patterns of disease progression, is an effective method to overcome current treatment bottlenecks and achieve personalized and precise intervention. Spatiotemporally specific drug targets refer to biomolecules whose expression or activity changes significantly under specific temporal or spatial conditions. These targets enable precise intervention in drug development. Such targets not only possess spatial specificity but also temporal specificity related to disease progression, typically nested within specific cell-cell interactions, metabolic states, immune signaling, or epigenetic regulatory contexts.

[0005] Although the research value of this field is widely recognized, computational methods capable of systematically reconstructing the characteristics of disease tissue microenvironments and restoring their spatiotemporal changes throughout the disease course remain relatively scarce. Existing models and algorithms suffer from two main problems: some spatial deconvolution methods heavily rely on high-cost technologies such as spatial transcriptomics, making them difficult to popularize; others, such as classic gene modularization algorithms like WCGNA and SCENIC, not only fail to integrate multi-dimensional reconstruction algorithms into the same computational biology framework but also fail to attribute the specific biological significance of individual modules.

[0006] Single-cell transcriptome sequencing (scRNA-seq) is a well-developed and widely used omics technology, with a large amount of single-cell sequencing data accumulated for metabolic chronic diseases. Single-cell transcriptome sequencing constructs the gene expression profile of each cell at the individual cell level, aiming to reveal gene expression levels within a single cell and understand cellular heterogeneity and functional diversity. However, this sequencing technology requires enzymatic digestion of the studied tissue, resulting in the loss of spatial location information of cells in the suspension. This spatial information is crucial for understanding tissue structural heterogeneity, cellular components, and interactions between neighboring cells within the tissue. Therefore, how to reconstruct the disease tissue microenvironment, elucidate disease mechanisms, and discover targets using single-cell transcriptome sequencing data has become a core technological challenge that urgently needs to be overcome in this field. Summary of the Invention

[0007] To address the technical problems existing in the prior art, embodiments of the present invention provide a method, system, and computer device for screening drug targets with spatiotemporal specificity. The technical solution is as follows:

[0008] A method for screening drug targets with spatiotemporal specificity, the method being based on single-cell transcriptome sequencing data, the method comprising:

[0009] (1) Quantitatively reconstructing the spatial location and functional patterns of cells in tissues, including:

[0010] 1.1) The single-cell transcriptome sequencing data is preprocessed, and the preprocessing includes the following steps:

[0011] 1.1.1) Data quality control and integration: The single-cell transcriptome sequencing data were subjected to quality control and normalization to remove low-quality cells and low-expression genes, and the data were integrated, principal component analysis was used for dimensionality reduction, and Louvain unsupervised clustering was performed.

[0012] 1.1.2) Cell type annotation: The unified manifold approximation projection algorithm is used to perform nonlinear dimensionality reduction on the clustering results and to annotate the cell types.

[0013] 1.2) Reconstruction of the spatial localization of a single cell, which includes the following steps:

[0014] 1.2.1) Construction of spatially localized gene sets: Based on spatial omics data, the tissue under study is divided into multiple regions, and gene sets with the highest expression in each region are constructed.

[0015] 1.2.2) Non-negative matrix factorization: The gene expression matrix of the total cell population obtained in step 1.1) is decomposed into the product of the marker gene set expression matrix corresponding to each spatial partition and the distribution score matrix of cells in different spatial partitions using the non-negative matrix factorization (NMF) Brunet algorithm.

[0016] 1.2.3) Perform UMAP dimensionality reduction and visualization on the matrix extracted from NMF to reconstruct the spatial localization of each cell subpopulation;

[0017] 1.3) Reconstruction of single-cell biological functional patterns, which includes the following steps:

[0018] 1.3.1) Single-sample gene set enrichment analysis ssGSEA: For each pathway in the specified functional module, ssGSEA enrichment scoring is performed based on the gene expression matrix of the total cell population obtained in step 1.1);

[0019] 1.3.2) Perform UMAP nonlinear dimensionality reduction, K-nearest neighbor graph construction, and Leiden unsupervised clustering on the ssGSEA score matrix;

[0020] 1.3.3) Output the classification results of specific functional characteristics of each single cell, and the enrichment scores of each functional pattern in each pathway;

[0021] (2) Screening for drug targets with spatiotemporal specificity includes the following steps:

[0022] 2.1) Perform cell-cell communication analysis, which includes: analyzing the communication probability and ligand-receptor interaction pairs between each cell subpopulation annotated in step 1.1) and other cell types, and establishing a mapping relationship between spatial / functional characteristics and cell interaction networks for the spatial partitions and functional patterns reconstructed for each cell in steps 1.2) and 1.3);

[0023] 2.2) Construct a gene regulatory network (GRN) centered on the microenvironment state of a specific tissue, which includes: calculating the interaction records between cell subpopulations and performing gene regulatory network analysis on hub genes;

[0024] 2.3) Target discovery, which includes: screening key regulatory factors and signaling nodes from the gene regulatory network GRN, and the obtained hub genes and their connected core nodes are potential spatiotemporally specific drug targets.

[0025] Optionally, in step 1.1.1), the quality control and normalization steps are as follows: use Cell Ranger to perform quality control and re-attachment of single-cell transcriptome sequencing data, perform Seurat quality control process on the output feature counting matrix, and use the DoubletFinder package to identify and remove cell duplexes;

[0026] And / or, in step 1.1.1), the steps of data integration, principal component analysis dimensionality reduction, and Louvain unsupervised clustering are as follows: IntegrateData is called to integrate the single-cell transcriptome sequencing data of different samples that have completed quality control and normalization; ScaleData and RunPCA are called to perform linear dimensionality reduction on the integrated Seurat object; the FindNeighbors function is called to construct the K nearest neighbor map to establish a cell similarity network; and the FindClusters function is called to perform unsupervised clustering using the Louvain algorithm.

[0027] And / or, in step 1.1.2), the steps for performing nonlinear dimensionality reduction are as follows: the clustering results are nonlinearly reduced using the UMAP algorithm, and the clustering results and grouping information are visualized using the CellDimPlot function in the SCP package;

[0028] And / or, in step 1.1.2), the steps for cell type annotation are as follows: refer to the cell marker database, complete the cell type annotation for all clusters, and use the FeaturePlot function to examine the expression distribution of marker genes in each cluster.

[0029] Optionally, in step 1.2.1), the structural characteristics of the tissue or organ under study are identified and it is divided into multiple partitions; the spatial omics data are processed to identify the gene expression status of each partition and the partition where each gene shows a peak expression level, and the R package or R function is used to organize them into a set of genes that show a peak expression level in each partition;

[0030] And / or, in step 1.2.2), the NMF package is used to perform non-negative matrix factorization on the total cell population gene expression matrix obtained in step 1.1). If the number of partitions is n, then the factorization rank is set to rank=n+1.

[0031] And / or, in step 1.2.3), the CreateDimReducObject function is used to embed the score matrix of each cell in the principal component direction into the Seurat object corresponding to the gene expression matrix of the total cell population. FindNeighbors and FindClusters are used in sequence to perform UMAP nonlinear dimensionality reduction and Louvain clustering on the NMF output again. The CellDimPlot function is used to display the spatial information reconstruction result in the two-dimensional feature space of UMAP based on the NMF dimensionality reduction result.

[0032] Optionally, in step 1.3.1), the ssGSEA enrichment scoring steps are as follows: Retrieve all functional pathways of the studied species and filter relevant pathway information; obtain the Entrez IDs of genes contained in each pathway and construct a gene set list for the studied pathways; use the ssGSEA algorithm in the GSVA package to perform gene set enrichment analysis on the gene expression matrix of the total cell population, and each cell will obtain an enrichment score for each gene set: ES i (S), and finally obtain the activity score matrix of each cell on the studied pathway;

[0033] And / or, in step 1.3.2), the activity score matrix of each cell on the studied pathway is reduced to UMAP dimension using the umap function in the uwot package, a KNN network is constructed based on the UMAP feature space using the FNN package, and unsupervised clustering is performed on the graph structure using the Leiden algorithm;

[0034] And / or, in step 1.3.2), it also includes: summarizing the clustering results by functional mode, calculating the average score of each functional mode on each pathway, and using the pheatmap package to draw a heatmap to show the activity score distribution of different functional modes;

[0035] And / or, in step 1.3.3), the UMAP dimensionality reduction result object is converted into a data frame using the as.data.frame function, and the UMAP dimensionality reduction result is visualized in a two-dimensional feature space using ggplot2 to create a scatter plot; the activity score matrix of a single cell on the studied pathway is organized by cluster using the group_by function, and the summarise function is used to calculate the mean ES value of each pathway for each cluster;

[0036] And / or, in step 1.3.3), it also includes: using the recode function to integrate the functional feature classification results of each single cell into the Seurat object obtained in step 1.1) or step 1.2).

[0037] Optionally, the cell-cell communication analysis is based on the CellChat package;

[0038] And / or, the steps for analysis using the CellChat package are as follows: construct a CellChat object from the Seurat object of the total cell population constructed in step 1.1), use the database provided by CellChat to annotate communication pathways, and sequentially perform gene expression enrichment analysis, ligand-receptor pair screening, communication probability estimation, pathway-level communication score calculation, and network aggregation.

[0039] Optionally, the construction of the gene regulatory network with a specific subpopulation as the core uses the CellChat package to calculate the interaction records between cell subpopulations and the RcisTarget package to perform gene regulatory network analysis on the hub gene.

[0040] And / or, the steps for performing gene regulatory network analysis on hub genes using the RcisTarget package are as follows: Load the motif ranking database and motif-TF annotation table of the species to which the CellChat database belongs, and perform motif enrichment calculation on the hub genes; screen for enriched motifs with NES scores ≥3, and extract high-confidence transcription factors and their regulated target genes; after standardization, use the ggraph package to draw the TF-Target network relationship diagram between transcription factors and target genes.

[0041] Optionally, the steps for screening key regulatory factors and signaling nodes from GRNs are as follows: The filter function is used to screen ligand-receptor pairs in the communication matrix constructed by CellChat that satisfy p < 0.05, and the communication path names, probabilities, and significance are extracted; the summary function is used to calculate the frequency and average communication strength of each ligand-receptor pair, and the two are combined into a candidate gene pair score matrix; a screening threshold is set for frequencies above the top 30 percentile and average communication strength above 0.01, and the screening results are the hub gene list. The key transcription factors obtained from RcisTarget analysis are then combined to identify the key regulatory factors and signaling nodes.

[0042] And / or, input the key regulatory factors and signal nodes into the drug-target database for retrieval, export the retrieval results, and visualize the undirected drug-target network using the igraph package; based on the drugs corresponding to the known targets included in the database, if the current indication of the drug is different from the disease being studied, then the repurposing of this drug can be considered in the future; if no known drug has been found for this target, then new drug development for this target can be considered in the future.

[0043] A screening system for drug targets with spatiotemporal specificity, the screening system being used to implement the aforementioned screening method, the screening system comprising:

[0044] (1) Quantitative reconstruction of the spatial localization and functional pattern of cells in tissues, including:

[0045] The data preprocessing module is used to preprocess the single-cell transcriptome sequencing data. The data preprocessing includes: quality control and normalization of the single-cell transcriptome sequencing data, data integration, principal component analysis dimensionality reduction and Louvain unsupervised clustering, and cell type annotation.

[0046] The spatial positioning reconstruction module is used to reconstruct the spatial positioning of single cells. The spatial positioning reconstruction includes: constructing a spatial positioning gene set, non-negative matrix factorization, and reconstructing the spatial positioning of each cell subpopulation.

[0047] The functional phenotype reconstruction module is used to reconstruct the functional patterns of single-cell biology. The functional phenotype reconstruction includes: single-sample gene set enrichment analysis, UMAP nonlinear dimensionality reduction, K nearest neighbor map construction and Leiden unsupervised clustering, and output of functional feature classification results and enrichment scores.

[0048] (2) A spatiotemporally specific drug target discovery module, which includes:

[0049] Cell-to-cell communication analysis module, used for performing cell-to-cell communication analysis;

[0050] A module for constructing gene regulatory networks, used to build GRN gene regulatory networks centered on specific tissue microenvironment states; and,

[0051] The target discovery module is used to screen key regulatory factors and signaling nodes from the gene regulatory network GRN. The obtained hub genes and their connected core nodes are potential spatiotemporally specific drug targets.

[0052] A computer device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the filtering method.

[0053] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the filtering method.

[0054] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following:

[0055] 1. Reduce the cost of omics experiments

[0056] Some spatial deconvolution methods heavily rely on high-cost technologies such as spatial transcriptomics, which are difficult to popularize due to limitations in experimental techniques and equipment. In contrast, single-cell RNA sequencing technology is relatively mature and less expensive. This invention can reconstruct spatial and biological functional feature patterns using only scRNA-seq data and known spatial and functional characteristics.

[0057] 2. The analysis results can be traced back to biological significance from multiple dimensions.

[0058] Classical gene modularization algorithms such as WCGNA and SCENIC have significant limitations in terms of data. They cannot integrate multi-dimensional reconstruction algorithms into the same computational biology framework, nor can they attribute the specific biological significance of a single module. This invention can characterize cell gene expression patterns from the perspectives of spatial distribution and multiple biological functional modes, and can also trace the biological significance of different patterns based on their characteristics, providing better guidance for clinical research in biology and basic medicine.

[0059] 3. Supports custom biological function analysis and non-model organism analysis.

[0060] Conventional spatial deconvolution and gene modularization algorithms often fail to selectively observe custom functional features, while this invention enables quantitative gene expression pattern characterization of arbitrary partitions and functional pathways of interest to the user; and compared to the high dependence of conventional methods on specific model organism databases, this invention supports the analysis of non-model organisms.

[0061] 4. Supports the screening of key target genes by combining drug-target databases, providing new ideas for drug repurposing or new drug development.

[0062] After constructing a gene regulatory network using this invention, the hub gene can be compared with existing drug-target databases to screen for therapeutic targets with the potential for drug repurposing or new drug development.

[0063] 5. Provide new research tools for basic research.

[0064] This invention can serve as a pre-analysis or auxiliary analysis tool in high-throughput omics research, and it complements other omics technologies (such as metabolomics, spatial omics, and epigenomics). On the one hand, this invention can be used to preliminarily identify potential functional hotspots at the single-cell transcriptional level, helping users narrow down their target scope, optimize sample design, and reduce experimental costs before spatial sequencing or metabolomics sequencing. On the other hand, its prediction results can also be cross-validated with actual metabolic throughput and spatial expression maps, thereby enhancing the credibility and interpretive depth of biological discoveries.

[0065] In summary, the technical solution of this invention has significant beneficial effects in reducing the cost of omics experiments and tracing the biological significance from multiple dimensions. Furthermore, this invention has broad applicability, supporting basic researchers in analyzing specific biological functions and non-model organisms, and supporting docking with drug-target databases. The method of this invention not only proposes and preliminarily verifies an innovative single-cell data analysis method, but also provides new ideas and technical pathways for drug development in metabolic diseases and other systemic diseases, possessing the foundation for further promotion and application. Attached Figure Description

[0066] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0067] Figure 1 This is a flowchart of the single-cell spatial localization and biological function reconstruction module provided in this embodiment of the invention;

[0068] Figure 2 This is a flowchart of the spatiotemporally specific drug target discovery module provided in this embodiment of the invention.

[0069] Figure 3 This is a diagram showing the arrangement of cells in a feature space after spatial localization reconstruction provided by an embodiment of the present invention, demonstrating the restoration of the spatial distribution of cells by the calculation method;

[0070] Figures 4A to 4C This is a diagram showing the arrangement of cells in the feature space after the reconstruction of the glucose and lipid metabolism functional pattern provided in this embodiment of the invention, demonstrating the restoration of the cellular glucose and lipid metabolism functional pattern by POALRIS. Figure 4A This is a UMAP dimensionality reduction diagram of the five metabolic modes; Figure 4B This is a graph showing the average ssGSEA scores of the five metabolic modes in each pathway. Figure 4C This is a reduced-dimensionality map of hepatocytes using UMAP (colored according to metabolic pattern, cell subpopulation, and modeling time).

[0071] Figures 5A to 5C This is a diagram showing the arrangement of cells in the feature space after the reconstruction of epigenetic gene expression patterns provided in this embodiment of the invention, demonstrating the restoration of cellular epigenetic gene expression patterns by POALRIS. Figure 5A This is a graph showing the average ssGSEA scores of the five epigenetic patterns in each pathway. Figure 5B This is a UMAP dimensionality reduction diagram of the five epigenetic patterns; Figure 5CThis is a reduced-dimensionality plot of hepatocytes using UMAP (colored according to epigenetic pattern, cell subset, and modeling time).

[0072] Figures 6A to 6D This is a diagram showing the overall immune interaction strength of cells grouped in different tissue microenvironments according to an embodiment of the present invention. Figure 6A This diagram shows how hepatocyte subsets are classified into five categories based on microenvironment information. Figure 6B It is a graph showing the total intensity of immune interactions involving hepatocyte subsets in different microenvironment groups; Figure 6C It is a graph showing the total intensity of immune interactions in different liver lobule regions during the statistical modeling process; Figure 6D This is a graph showing the relationship between immune signal emission and reception and the signal intensity among different cell types;

[0073] Figures 7A to 7C This is a diagram illustrating the bias in the emission / reception of inflammation-related molecular signals by cells under different tissue microenvironment groupings provided in this embodiment of the invention. Figure 7A This is a diagram showing the central changes in pro-inflammatory / regulatory immune interactions during the modeling process; Figure 7B This is a diagram showing the pro-inflammatory / regulatory bias in the interaction between dendritic cells and hepatocytes of different microenvironment classifications; Figure 7C It is Mlxipl + A graph showing the statistical bias of pro-inflammatory / regulatory tendencies in hepatocytes;

[0074] Figures 8A to 8C This is a diagram showing the alignment of gene regulatory networks and potential targets in key molecular processes provided in this invention with data from a database. Figures 8A to 8B It is dendritic cells, Mlxipl + A diagram of GRN, which is constructed around hepatocyte-mediated immune interactions; Figure 8C Based on Mlxipl + Key genes involved in the process of hepatocytes sending pro-inflammatory signals to dendritic cells are identified in a reported drug-target map retrieved from DGIdb. Detailed Implementation

[0075] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0076] Definitions of abbreviations, their English equivalents, and key terms:

[0077] POLARIS (Pathway-Oriented Layered Assembly for Regulatory InferenceSystems);

[0078] scRNA-seq: Single-cell transcriptome sequencing;

[0079] ssGSEA: Single-sample gene set enrichment analysis.

[0080] To address the challenges of existing technologies, this invention proposes a method for multidimensional reconstruction of the microenvironment and drug target identification based on scRNA-seq data. This method enables quantitative reconstruction of the spatial localization and functional patterns of cells within tissues, and the discovery of spatiotemporally specific drug targets. It offers advantages such as applicability to non-model species, support for drug repurposing and new target development, and provides new tools and approaches for mechanistic research and targeted therapy of complex systemic diseases such as metabolic diseases and tumors.

[0081] This invention first provides a method for screening spatiotemporally specific drug targets based on single-cell transcriptome data. The method includes quantitatively reconstructing the spatial localization and functional patterns of cells in tissues, and discovering spatiotemporally specific drug targets.

[0082] The quantitative reconstruction of the spatial localization and functional patterns of cells in tissues includes (1) data preprocessing; (2) reconstruction of the spatial localization of single cells and (3) reconstruction of the biological functional patterns of single cells.

[0083] The data preprocessing in (1) includes:

[0084] (1.1) Data quality control and integration: Quality control and normalization of single-cell transcriptome sequencing data of biological samples, removal of low-quality cells and low-expression genes, and data integration, principal component analysis (PCA) dimensionality reduction and Louvain unsupervised clustering for multiple samples;

[0085] (1.2) Cell type annotation: The clustering results were nonlinearly reduced using the Unified Manifold Approximation Projection (UMAP) algorithm, and cell type annotation was performed by consulting literature and databases.

[0086] The reconstruction of single-cell spatial localization mentioned in (2) includes:

[0087] (2.1) Construction of spatially localized gene sets: Based on the literature on spatial transcriptomics or known spatial omics data, the tissue under study was divided into multiple regions, and gene sets with the highest expression in each region were constructed respectively;

[0088] (2.2) Non-negative matrix factorization: Extract the preprocessed scRNA-seq data, and use the non-negative matrix factorization (NMF) Brunet algorithm to decompose the gene expression matrix of the total cell population into the product of the gene set matrix of different partitions and the distribution weight matrix of individual cells in different spatial partitions;

[0089] (2.3) Perform UMAP dimensionality reduction and visualization on the matrix extracted from NMF to reconstruct the spatial localization of each cell subpopulation.

[0090] The reconstruction of the single-cell biological functional pattern mentioned in (3) includes:

[0091] (3.1) Single-sample gene set (ssGSEA) enrichment analysis: For each pathway in the specified functional module, ssGSEA enrichment scoring is performed based on the single-cell expression matrix;

[0092] (3.2) Perform UMAP nonlinear dimensionality reduction on the ssGSEA score matrix, construct the K nearest neighbor graph (KNN graph) and perform Leiden unsupervised clustering;

[0093] (3.3) Output the classification results of specific functional characteristics of each single cell, and the enrichment scores of each functional pattern in each pathway.

[0094] The spatiotemporally specific drug target discovery method includes (1) cell-cell communication analysis; (2) constructing a gene regulatory network (GRN) with a specific subgroup as the core; and (3) target discovery.

[0095] In some embodiments, the cell-cell communication analysis is based on the CellChat package, which analyzes the communication probability and ligand-receptor interaction pairs between each cell subpopulation and other cell types at different spatiotemporal scales, thereby establishing a mapping relationship between the spatial / functional features characterized by this method and the cell interaction network.

[0096] In some embodiments, the construction of a gene regulatory network centered on a specific subpopulation is performed by calculating interaction records between cell subpopulations using the CellChat package workflow, and the RcisTarget package is used to perform gene regulatory network analysis on the hub gene.

[0097] In some embodiments, target discovery involves screening key regulatory factors and signaling nodes from GRNs. These hub genes and their connected core nodes represent potential spatiotemporally specific drug targets, providing a theoretical basis and experimental evidence for subsequent drug interventions.

[0098] The present invention also provides applications of the above method, wherein the applications are a) to perform analysis of specific tissue microenvironments of systemic diseases such as metabolic diseases and tumors or other biological samples; or b) to discover spatiotemporally specific drug targets in the above samples; or c) to study key targets in cell differentiation and individual development processes, wherein the applications are for non-diagnostic or non-therapeutic purposes.

[0099] The present invention also provides a script-integrated computation toolkit, which has the function of implementing the computation method.

[0100] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0101] Example 1

[0102] A method for screening drug targets with spatiotemporal specificity, comprising:

[0103] 1. A method for quantitative reconstruction of spatial localization and functional patterns of cells in tissues based on single-cell transcriptome data. This method can reconstruct the spatial localization of single cells in tissues and gene expression patterns of specific biological functions from scRNA-seq data. The overall workflow framework is as follows: Figure 1 As shown. The specific technical solution includes the following steps:

[0104] 1.1 Data Preprocessing:

[0105] (1.1.1) Data quality control and integration: The single-cell transcriptome sequencing data of biological samples were subjected to quality control and normalization to remove low-quality cells and low-expression genes. For scRNA-seq data of samples with different experimental treatments, after quality control and normalization, data integration, principal component analysis (PCA) dimensionality reduction and Louvain unsupervised clustering were performed.

[0106] In one optional embodiment, the specific method for quality control and normalization is as follows: Cell Ranger (v 9.0.0) is used for quality control and backtracking of raw scRNA-seq reads. The output feature count matrix is ​​then subjected to a Seurat quality control workflow (Seurat 5.1.0), with the following selection criteria: the number of genes detected per cell (nFeature_RNA) is between 200 and 5000, and the total number of transcripts (nCount_RNA) is less than 20000. The Seurat functions NormalizeData, FindVariableFeatures, ScaleData, and RunPCA are used to standardize the data and select features, choosing the top 30 principal components. Subsequently, the DoubletFinder package (v 2.0.4) is used to identify and remove possible cell doublets. The expected number of doublets is simulated based on the 10x Genomics recommended formula (setting a double-cell rate of 0.8% per thousand cells), retaining only cells identified as single cells for subsequent analysis. The advantage of this step is that it provides a more comprehensive cleaning of the raw data, removing low-quality cells and duplexes, thus ensuring the quality of data for subsequent analysis.

[0107] In one alternative embodiment, the specific methods for integrating sample data from different experimental treatments, principal component analysis (PCA) dimensionality reduction, and Louvain unsupervised clustering are as follows:

[0108] For scRNA-seq data of samples from different experimental treatments that have undergone quality control and normalization, FindIntegrationAnchors was used to identify anchor points. The top 30 principal components were selected as integration dimensions, and the top 2000 hypervariable genes were used as anchor features. IntegrateData was called to integrate multiple samples, resulting in a new Seurat object. This Seurat object was then subjected to linear dimensionality reduction using ScaleData and RunPCA. Finally, the FindNeighbors function was called to extract the top n dimensions to construct a K-nearest neighbor graph (KNN), thereby establishing a cell similarity network (the ElbowPlot function can be used to draw an elbow plot on the Seurat object, with the optimal value of n corresponding to the corner of the image). Based on this, the FindClusters function was called with a resolution of 0.8, and the Louvain algorithm was used for unsupervised clustering. After completing the above integration, dimensionality reduction, and clustering steps, the gene expression matrix of this Seurat object is called the total cell population gene expression matrix. This matrix contains the gene expression information of all cells obtained under different experimental treatments in this study.

[0109] The advantage of performing multi-group sample data integration, principal component analysis (PCA) dimensionality reduction, and Louvain unsupervised clustering steps is that it integrates scRNA-seq data of biological samples from different experimental treatments into the same matrix, which facilitates subsequent unified analysis.

[0110] (1.1.2) Cell type annotation: The clustering results were nonlinearly reduced using the Unified Manifold Approximation Projection (UMAP) algorithm, and cell type annotation was performed by consulting literature and databases.

[0111] In one alternative embodiment, the steps for nonlinear dimensionality reduction are as follows: the clustering results are nonlinearly reduced using the UMAP algorithm (selecting the first 30 principal components), and both the clustering results and grouping information are visualized using the CellDimPlot function in the SCP package (v 0.5.6).

[0112] In one optional embodiment, the cell type annotation steps are as follows: Cell type annotation for all clusters is completed by comprehensively reviewing and referencing cell marker databases and literature. Simultaneously, the FeaturePlot function is used to examine the expression distribution of marker genes in each cluster. The final annotation results are manually mapped and assigned to the CellType field of meta.data, and this field is used as the primary identifier for subsequent analysis and visualization.

[0113] The advantage of performing cell type annotation is that it maps the gene expression of cells to the corresponding cell types, which facilitates further research based on the biological functions of the corresponding cell types.

[0114] 1.2 Reconstruction of the spatial localization of a single cell:

[0115] (1.2.1) Construction of spatially localized gene sets: Based on the literature on spatial transcriptomics or known spatial omics data, the tissue under study is divided into multiple regions, and gene sets with the highest expression in each region are constructed respectively;

[0116] In an optional embodiment, the steps of dividing the tissue into multiple partitions and constructing the gene set with the highest expression in each region are as follows: The structural characteristics of the tissue or organ under study are clarified by comprehensively reviewing databases and literature, and based on this prior knowledge, it is divided into multiple partitions (e.g., liver lobules are divided into portal vein region, intermediate region, and central vein region); simultaneously, the spatial omics data used as spatial partition marker genes are processed to clarify the gene expression status of each partition, as well as the partition where each gene shows a peak expression level, and the gene set with the peak expression level in each partition is organized using R packages such as dplyr and R functions such as subset, which.max, and filter.

[0117] (1.2.2) Nonnegative matrix factorization: The gene expression matrix of the total cell population obtained in step (1.1) (e.g., the gene expression matrix of the total cell population obtained in step (1.1.2)) is decomposed into the product of the marker gene set expression matrix corresponding to each spatial partition and the distribution score matrix of cells in different spatial partitions using the nonnegative matrix factorization (NMF) Brunet algorithm.

[0118] Nonnegative Matrix Factorization (NMF) is a matrix factorization method proposed by Lee and Seung in *Nature* in 1999. It requires all components to be nonnegative and achieves non-linear dimensionality reduction. This method captures local features of data through pure additive description, its psychological basis stemming from the property that overall perception is composed of partial perceptions. Sparse representation enhances the rationality of data interpretation. From a multivariate statistical perspective, NMF simplifies high-dimensional data into low-dimensional patterns while preserving information; from an algebraic perspective, it reveals the inherent nonnegative decomposition form of the data. This method is widely used in image analysis, text clustering, speech processing, and other fields, and its nonnegativity constraint effectively adapts to the practical needs of physical signal and data processing.

[0119] In an optional embodiment, the step of decomposing the product of the marker gene set expression matrix corresponding to each spatial partition and the cell distribution score matrix in different spatial partitions is as follows: Let the gene expression matrix of the total cell population be... Where X represents the gene expression matrix of the total cell population, Let X be an m-row, n-column matrix, where m represents the number of cells (encoded by cell barcodes in the actual matrix), and n represents the total number of marker genes in different partitions. NMF decomposes this matrix into the product of two non-negative matrices:

[0120]

[0121] in, Let W represent the distribution score matrix of cells across different spatial partitions. Let W be an m x k ​​matrix, where m represents the number of cells (encoded by cell barcodes in the actual matrix), and k represents the number of spatial partitions. The expression of marker gene sets corresponding to each spatial partition is represented by H, where H represents the marker gene set expression matrix corresponding to each spatial partition. H represents a k-row, n-column matrix, where n represents the total number of marker genes in different partitions, and k represents the number of spatial partitions.

[0122] In one optional embodiment, the NMF package (v 0.24.0) is used to perform non-negative matrix factorization on the gene expression matrix of the total cell population. If the number of partitions is n, the factorization rank is set to rank=n+1. The Brunet algorithm is used to obtain the score matrix of each cell in the principal component direction of the space (this score matrix is ​​the transpose of the coefficient matrix).

[0123] The advantage of performing nonnegative matrix factorization is that, in principle, nonnegative matrix factorization belongs to semi-supervised dimensionality reduction algorithms. Compared with traditional unsupervised dimensionality reduction algorithms, it has higher result interpretability and a larger user-defined space. It can selectively extract low-dimensional information from the high-dimensional matrix: the gene expression matrix of the total cell population, from the user-defined spatial partition marker genes: the cell's localization tendency in different partitions.

[0124] (1.2.3) Perform UMAP dimensionality reduction and visualization on the matrix extracted from NMF to reconstruct the spatial localization of each cell subpopulation.

[0125] In an optional embodiment, the steps for UMAP dimensionality reduction, visualization, and reconstruction of the spatial localization of each cell subpopulation are as follows: The `CreateDimReducObject` function is used to embed the score matrix of each cell along the principal component direction into a Seurat object corresponding to the gene expression matrix of the total cell population. Then, `FindNeighbors` (the number of dimes should be consistent with the linear dimensionality reduction in (1.1.1)) and `FindClusters` (different numbers of clusters can be obtained by adjusting the resolution) are used sequentially to perform UMAP nonlinear dimensionality reduction and Louvain clustering on the NMF output. Finally, the `CellDimPlot` function is used to display the spatial information reconstruction results in the two-dimensional feature space of UMAP, based on the NMF dimensionality reduction results.

[0126] 1.3. Reconstruction of Functional Patterns in Single-Cell Biology:

[0127] (1.3.1) Single Sample Gene Set (ssGSEA) Enrichment Analysis: For each pathway in the specified functional module, ssGSEA enrichment scoring is performed based on the gene expression matrix of the total cell population obtained in step 1.1.

[0128] In an optional embodiment, the ssGSEA enrichment scoring steps are as follows: All functional pathways of the studied species are retrieved from the functional pathway database, and relevant pathway information is filtered. The Entrez IDs of genes contained in each pathway are obtained using the KEGGREST package (v 1.42.0), and converted to gene symbol form using the bioomaRt package to construct a gene set list for the studied pathway. Subsequently, the ssGSEA algorithm in the GSVA package (v 1.50.5) is used to perform gene set enrichment analysis on the gene expression matrix of the total cell population. Each cell receives an enrichment score (ES) for each gene set: ES i (S), ultimately yielding an activity score matrix for each cell on the studied pathway.

[0129] The advantages of performing single-sample gene set (ssGSEA) enrichment analysis are twofold: firstly, in terms of sample quantity, ssGSEA is suitable for enriching a large number of samples by treating a single cell as a single sample; secondly, in terms of data type, compared to the GSVA enrichment algorithm, ssGSEA is more stable for sparse matrices such as scRNA-seq gene expression matrices.

[0130] (1.3.2) Perform UMAP nonlinear dimensionality reduction on the ssGSEA score matrix, construct the K nearest neighbor graph (KNN graph), and perform Leiden unsupervised clustering;

[0131] In an optional embodiment, the steps for performing UMAP nonlinear dimensionality reduction, K-nearest neighbor graph (KNN) construction, and Leiden unsupervised clustering are as follows: For each cell's activity score matrix on the studied pathway, UMAP dimensionality reduction is performed using the umap function in the uwot package (v 0.2.3) (setting n_neighbors = 30, min_dist = 0.3). Let U represent the low-dimensional ES matrix of gene sets with different functional modules after dimensionality reduction by UMAP. U represents an m x d matrix, where m is the number of cells (encoded by cell barcodes in the actual matrix) and d is the number of functional pathways studied. A KNN network (K = 50) is constructed based on the UMAP feature space using the FNN package (v 1.1.4.1). Unsupervised clustering is performed on the graph structure using the Leiden algorithm (the resolution parameter can be set according to the desired number of clusters in the actual application). The clustering results are summarized by functional mode, and the average score for each functional mode on each pathway is calculated. A heatmap is generated using the pheatmap package (v1.0.12) to display the activity score distribution of different functional modes.

[0132] The advantage of this step is that it allows for the visualization of the heterogeneity of all cells in terms of functional patterns through nonlinear dimensionality reduction. Furthermore, the Leiden clustering algorithm provides more reliable cluster partitioning compared to the conventional Louvain algorithm, thus avoiding fragmented functional clusters.

[0133] (1.3.3) Output the classification results of specific functional characteristics of each single cell, and the enrichment scores of each functional pattern in each pathway.

[0134] In one optional embodiment, the steps for outputting the classification results of specific functional features for each single cell are as follows: The UMAP dimensionality reduction result object is converted into a data frame using the `as.data.frame` function, and a scatter plot is created using ggplot2 to visualize the UMAP dimensionality reduction result in a two-dimensional feature space. The steps for outputting the enrichment scores of each functional pattern in each pathway are as follows: The activity score matrix of a single cell in the studied pathway is organized by cluster using the `group_by` function, and the mean ES value for each pathway is calculated using the `summarise` function for each cluster (i.e., each functional pattern).

[0135] The `recode` function is used to integrate the functional feature classification results of each single cell into the Seurat object obtained in step 1.1 or 1.2 (depending on whether the spatial localization information of the single cell is reconstructed) for further analysis.

[0136] 2. A spatiotemporally specific drug target discovery method. This method can perform cell-cell communication analysis on cell subpopulations located in specific tissue microenvironments as calculated above, and perform gene regulatory network analysis on hub genes in key molecular processes. Hub genes and core node molecules in the regulatory network are potential targets. The overall process framework of the method is as follows: Figure 2 As shown. The specific technology includes the following steps:

[0137] 2.1 Cell-to-cell communication analysis

[0138] Cell-cell communication analysis is based on the CellChat package. It analyzes the communication probability and ligand-receptor interaction pairs between each cell subpopulation annotated in step 1.1 and other cell types. Combined with the spatial partitions and functional patterns reconstructed for each cell in steps 1.2 and 1.3, a mapping relationship between the spatial / functional features characterized by this method and the cell interaction network can be established.

[0139] After steps (1.2.3) and (1.3.3) are completed, the spatial information and functional pattern information have been embedded into the Seurat object corresponding to the gene expression matrix of the total cell population. At this point, it is possible to characterize the spatial distribution and functional pattern characteristics of each cell subpopulation annotated in step (1.1.2). The cell interaction network is calculated on this subpopulation as a unit. By organizing the cell communication analysis results through data processing R packages such as dplyr, the functional pattern and spatial partitioning of each cell subpopulation as well as ligand-receptor interaction information can be obtained. This information is the mapping relationship mentioned above.

[0140] Cell communication refers to the ability of cells to receive, process, and transmit signals to their environment and themselves; it is a fundamental property of all cells in every organism, such as bacteria, plants, and animals. Cell-to-cell communication mediated by ligand-receptor complexes is crucial for coordinating various biological processes, such as development, differentiation, and inflammation.

[0141] CellChat is an R package that can quantitatively infer and analyze intercellular communication networks from single-cell transcriptome sequencing (scRNA-seq) data. It requires cellular gene expression data as input and establishes the probability of cell-cell communication by integrating prior knowledge of the interactions between gene expression and signal ligands, receptors and their cofactors, thereby predicting intercellular communication and providing a variety of visualization methods.

[0142] In an optional embodiment, the steps for analysis using the CellChat package are as follows: A CellChat object (CellChat 2.1.2) is constructed from the Seurat object of the total cell population built in step 1.1, and communication pathway annotation is performed using the database provided by CellChat. If the species under study is human or mouse, the built-in CellChat database can be used directly; if the species under study is a non-model species, only homologous genes from closely related species can be retained for further analysis. Gene expression enrichment analysis, ligand-receptor pair screening, communication probability estimation, pathway-level communication score calculation, and network aggregation are performed sequentially, with all parameters set to default.

[0143] The advantage of this step compared to existing technologies is that it no longer performs cell communication analysis based solely on the single-cell transcriptome itself. Instead, it changes the perspective and analyzes individual cells within the spatial and functional characteristics reconstructed in steps 1.2 and 1.3, which can provide more biological significance for intercellular communication.

[0144] 2.2 Constructing a gene regulatory network (GRN) centered on a specific subgroup

[0145] A gene regulatory network centered on a specific subpopulation was constructed. Interaction records between cell subpopulations were calculated using the CellChat package, and the RcisTarget package was used to perform gene regulatory network analysis on the hub gene.

[0146] RcisTarget is an R package for constructing gene regulatory networks and analyzing transcription factors. It can identify potential transcription factor (TF) regulatory networks from a set of genes and infer the mechanisms of action of transcription factors through motif analysis. This package is primarily used in transcriptional regulation research, particularly for predicting upstream regulators of gene co-expression networks or differentially expressed gene sets.

[0147] In an optional embodiment, the steps for performing gene regulatory network analysis on hub genes using the RcisTarget package are as follows: Load the motif ranking database and motif-TF annotation table for the species to which the CellChat database belongs, and perform motif enrichment calculations on the hub genes. Screen for enriched motifs with NES scores ≥3, and extract high-confidence transcription factors (TF_highConf) and their regulated target genes. After standardization, use the ggraph package (v 2.2.1) to plot the TF-Target network relationship between transcription factors and target genes.

[0148] 2.3 Target Discovery

[0149] Key regulatory factors and signaling nodes are screened from GRNs. These hub genes and their connected core nodes are potential spatiotemporally specific drug targets, providing a theoretical basis and experimental evidence for subsequent drug interventions.

[0150] In an optional embodiment, the steps for screening key regulatory factors and signaling nodes from GRNs are as follows: Using the filter function, ligand-receptor pairs satisfying p < 0.05 in the communication matrix constructed by CellChat are screened as needed, and the communication pathway names, probabilities, and significance are extracted. Using the summary function, the frequency (i.e., the number of times it appears in different pathways or cell interactions) and average communication strength (mean prob value) of each ligand-receptor pair are calculated, and the two are combined into a candidate gene pair score matrix. A screening threshold is set at a frequency above the top 30 percentile and an average communication strength above 0.01. The screening results are the hub gene list (not limited to ligands or receptors). The key transcription factors obtained from RcisTarget analysis are then combined to identify the key regulatory factors and signaling nodes.

[0151] In one optional embodiment, the aforementioned key regulatory factors and signaling nodes are input into a drug-target database for retrieval. The retrieval results are exported, and the undirected drug-target network is visualized using the igraph package (v 2.0.3). Based on the drugs corresponding to known targets included in the database, if the current indication of the drug differs from the disease under investigation, the drug's "drug repurposing" can be considered: that is, targeting the same target and using the same drug to treat different diseases; if no known drug has been found for this target, new drug development for this target can be considered.

[0152] The advantage of this target discovery step compared to existing technologies is that it no longer relies solely on single-cell transcriptome analysis for cell communication analysis. Instead, it changes the perspective by placing individual cells within the spatial and functional features reconstructed in steps 1.2 and 1.3 for comprehensive analysis. This can provide more biological significance for the selected hub genes and signaling nodes.

[0153] Example 2

[0154] Practical application: The method of this invention can be applied to scRNA-seq data of biological samples after disease modeling to analyze corresponding spatiotemporally specific drug targets.

[0155] The input is the gene expression matrix of the total cell population, as well as the reported information of specific spatiotemporal and functional markers; the output is the spatial localization of each single cell and the classification results of specific functional gene expression patterns, as well as the enrichment scores of each functional pattern in each pathway and the cell-cell interaction and gene regulatory network of specific pattern cell subpopulations.

[0156] Validation by Example: Using scRNA-seq data from liver tissue of a Syrian hamster model of metabolic dysfunction-related steatohepatitis (MASH) at weeks 0, 2, and 6 (modeling was performed using a conventional / high-fat diet; liver tissue was harvested at weeks 0, 2, and 6 to prepare single-cell suspensions, and single-cell transcriptome sequencing was performed using an Illumina high-throughput sequencing platform), the method of this invention was applied for analysis. Using the spatial localization reconstruction module, a strong arrangement bias of annotated hepatocyte subsets in the liver lobule tissue space was observed in the UMAP-reduced two-dimensional feature space. Results are shown in [link to results]. Figure 3 Using the functional phenotype reconstruction module, strong heterogeneity was observed in the gene expression patterns of annotated hepatocyte subsets across different functions. See the results below. Figures 4A to 4C and Figures 5A to 5C Furthermore, both types of remodeling modules can be attributed to specific biological meanings. For example, hepatocytes in Mode 3 of glucose and lipid metabolism are in a state of highly active glucose and lipid metabolism; hepatocytes in Mode 5 of epigenetic-related gene expression are in a state of active chromatin remodeling, etc.

[0157] Target Discovery: The spatiotemporally specific drug target discovery module of this invention first performs cell-cell communication analysis based on annotated cell subpopulations. Based on a custom set of pro-inflammatory and anti-inflammatory molecules, the communication signal preferences of different tissue microenvironment cell subpopulations can be analyzed. Combined with the results of upstream analysis, it can be observed that the total intensity of immune interaction in the proximal portal vein region (Z1) of the liver lobule shows a trend of "first decreasing and then increasing"; the intermediate region (Z2) continuously increases over time, while the central venous region (Z3) remains at a low level. See [see details]. Figures 6A to 6D Dendritic cells (DCs) received the most pro-inflammatory signals from Group 2 (Abcb11⁺, Mlxipl⁺ hepatocytes), but only sent pro-inflammatory signals to Group 4 (Cd44⁺), suggesting a specific bias in immune interactions between different DC subsets. This analysis confirms the pro-inflammatory / anti-inflammatory functional bias of immune cells observed in basic experimental studies during disease development. Gene regulatory network analysis of the hub genes screened in the above process showed that the immune-related interactions of Mlxipl⁺ hepatocytes were relatively more complex. These hepatocytes are known to be mainly located in the region near the portal vein (Z1), a spatial localization characteristic consistent with the phenomenon of immune cell infiltration in the portal vein region during the course of MASH described in the literature. (See attached results). Figures 7A to 7C .

[0158] After characterizing the molecular behavior of this specific microenvironment, further target screening and drug repositioning analysis will be conducted to explore potential intervention methods specific to the tissue microenvironment. Regarding drug target identification, two scenarios are considered: first, there are already reported drugs for the target gene that have no known MASH applications, possessing the research potential for "drug repurposing"; second, there are no known corresponding drugs for the target gene, potentially representing entirely new therapeutic targets.

[0159] For the first scenario, the hub genes extracted from the aforementioned GRNs were input into the DGIdb database for retrieval to obtain information on genes covered by developed drugs and their corresponding drugs. A drug-target network map was then constructed, and the results are shown in [link to results]. Figures 8A to 8C The analysis results indicate that the drugs corresponding to these targets cover multiple therapeutic areas, including hormone replacement, anti-tumor, anticoagulation, and anti-diabetic applications.

[0160] In the second case, the number of edges of each hub gene in the constructed GRN was counted, and it was found that there are no known corresponding targeted drugs for the Rarres2 and Angptl4 genes, which have the potential to be novel targets.

[0161] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for screening drug targets with spatiotemporal specificity, characterized in that the screening method is based on single-cell transcriptome sequencing data, the method comprising: (1) Quantitatively reconstructing the spatial location and functional patterns of cells in tissues, including: 1.1) The single-cell transcriptome sequencing data is preprocessed, and the preprocessing includes the following steps: 1.1.1) Data quality control and integration: The single-cell transcriptome sequencing data were subjected to quality control and normalization to remove low-quality cells and low-expression genes, and the data were integrated, principal component analysis was used for dimensionality reduction, and Louvain unsupervised clustering was performed. 1.1.2) Cell type annotation: The unified manifold approximation projection algorithm is used to perform nonlinear dimensionality reduction on the clustering results and to annotate cell types. 1.2) Reconstruction of the spatial localization of a single cell, which includes the following steps: 1.2.1) Construction of spatially localized gene sets: Based on spatial omics data, the tissue under study is divided into multiple regions, and gene sets with the highest expression in each region are constructed. 1.2.2) Non-negative matrix factorization: The gene expression matrix of the total cell population obtained in step 1.1) is decomposed into the product of the marker gene set expression matrix corresponding to each spatial partition and the distribution score matrix of cells in different spatial partitions using the non-negative matrix factorization (NMF) Brunet algorithm. 1.2.3) Perform UMAP dimensionality reduction and visualization on the matrix extracted from NMF to reconstruct the spatial localization of each cell subpopulation; 1.3) Reconstruction of single-cell biological functional patterns, which includes the following steps: 1.3.1) Single-sample gene set enrichment analysis ssGSEA: For each pathway in the specified functional module, ssGSEA enrichment scoring is performed based on the gene expression matrix of the total cell population obtained in step 1.1); 1.3.2) Perform UMAP nonlinear dimensionality reduction, K-nearest neighbor graph construction, and Leiden unsupervised clustering on the ssGSEA score matrix; 1.3.3) Output the classification results of specific functional characteristics of each single cell, and the enrichment scores of each functional pattern in each pathway; (2) Screening for drug targets with spatiotemporal specificity includes the following steps: 2.1) Perform cell-cell communication analysis, which includes: analyzing the communication probability and ligand-receptor interaction pairs between each cell subpopulation annotated in step 1.1) and other cell types, and establishing a mapping relationship between spatial / functional characteristics and cell interaction networks for the spatial partitions and functional patterns reconstructed for each cell in steps 1.2) and 1.3); 2.2) Construct a gene regulatory network (GRN) centered on the microenvironment state of a specific tissue, which includes: calculating the interaction records between cell subpopulations and performing gene regulatory network analysis on hub genes; 2.3) Target discovery, which includes: screening key regulatory factors and signaling nodes from the gene regulatory network GRN, and the obtained hub genes and their connected core nodes are potential spatiotemporally specific drug targets.

2. The screening method according to claim 1, characterized in that, In step 1.1.1), the quality control and normalization steps are as follows: Cell Ranger is used to perform quality control and re-attachment of single-cell transcriptome sequencing data, Seurat quality control workflow is executed on the output feature counting matrix, and the DoubletFinder package is used to identify and remove cell duplexes; And / or, in step 1.1.1), the steps of data integration, principal component analysis dimensionality reduction, and Louvain unsupervised clustering are as follows: IntegrateData is called to integrate the single-cell transcriptome sequencing data of different samples that have completed quality control and normalization; ScaleData and RunPCA are called to perform linear dimensionality reduction on the integrated Seurat object; the FindNeighbors function is called to construct the K nearest neighbor map to establish a cell similarity network; and the FindClusters function is called to perform unsupervised clustering using the Louvain algorithm. And / or, in step 1.1.2), the steps for performing nonlinear dimensionality reduction are as follows: the clustering results are nonlinearly reduced using the UMAP algorithm, and the clustering results and grouping information are visualized using the CellDimPlot function in the SCP package; And / or, in step 1.1.2), the steps for cell type annotation are as follows: refer to the cell marker database, complete the cell type annotation for all clusters, and use the FeaturePlot function to examine the expression distribution of marker genes in each cluster.

3. The screening method according to claim 1, characterized in that, In step 1.2.1), the structural characteristics of the tissue or organ under study are identified and it is divided into multiple partitions; the spatial omics data are processed to identify the gene expression status of each partition and the partition where each gene shows a peak expression level. The R package or R function is used to organize the data into a set of genes that show a peak expression level in each partition. And / or, in step 1.2.2), the NMF package is used to perform non-negative matrix factorization on the total cell population gene expression matrix obtained in step 1.1). If the number of partitions is n, then the factorization rank is set to rank=n+1. And / or, in step 1.2.3), the CreateDimReducObject function is used to embed the score matrix of each cell in the principal component direction into the Seurat object corresponding to the gene expression matrix of the total cell population. FindNeighbors and FindClusters are used in sequence to perform UMAP nonlinear dimensionality reduction and Louvain clustering on the NMF output again. The CellDimPlot function is used to display the spatial information reconstruction result in the two-dimensional feature space of UMAP based on the NMF dimensionality reduction result.

4. The screening method according to claim 1, characterized in that, In step 1.3.1), the steps of ssGSEA enrichment scoring are as follows: retrieve all functional pathways of the studied species and filter relevant pathway information; obtain the Entrez ID of the genes contained in each pathway and construct a gene set list of the studied pathways; Using the ssGSEA algorithm from the GSVA package, gene set enrichment analysis is performed on the gene expression matrix of the total cell population. Each cell receives an enrichment score for each gene set: ES. i (S), and finally obtain the activity score matrix of each cell on the studied pathway; And / or, in step 1.3.2), the activity score matrix of each cell on the studied pathway is reduced to UMAP dimension using the umap function in the uwot package, a KNN network is constructed based on the UMAP feature space using the FNN package, and unsupervised clustering is performed on the graph structure using the Leiden algorithm; And / or, in step 1.3.2), it also includes: summarizing the clustering results by functional mode, calculating the average score of each functional mode on each pathway, and using the pheatmap package to draw a heatmap to show the activity score distribution of different functional modes; And / or, in step 1.3.3), the UMAP dimensionality reduction result object is converted into a data frame using the as.data.frame function, and the UMAP dimensionality reduction result is visualized in a two-dimensional feature space using ggplot2 to create a scatter plot; the activity score matrix of a single cell on the studied pathway is organized by cluster using the group_by function, and the summarise function is used to calculate the mean ES value of each pathway for each cluster; And / or, in step 1.3.3), it also includes: using the recode function to integrate the functional feature classification results of each single cell into the Seurat object obtained in step 1.1) or step 1.2).

5. The screening method according to claim 1, characterized in that, The cell-cell communication analysis was based on the CellChat package; And / or, the steps for analysis using the CellChat package are as follows: construct a CellChat object from the Seurat object of the total cell population constructed in step 1.1), use the database provided by CellChat to annotate communication pathways, and sequentially perform gene expression enrichment analysis, ligand-receptor pair screening, communication probability estimation, pathway-level communication score calculation, and network aggregation.

6. The screening method according to claim 1, characterized in that, The construction of a gene regulatory network centered on a specific subpopulation uses the CellChat package to calculate interaction records between cell subpopulations and the RcisTarget package to perform gene regulatory network analysis on the hub gene. And / or, the steps for performing gene regulatory network analysis on hub genes using the RcisTarget package are as follows: Load the motif ranking database and motif-TF annotation table of the species to which the CellChat database belongs, and perform motif enrichment calculation on the hub genes; screen for enriched motifs with NES scores ≥3, and extract high-confidence transcription factors and their regulated target genes; after standardization, use the ggraph package to draw the TF-Target network relationship diagram between transcription factors and target genes.

7. The screening method according to claim 1, characterized in that, The steps for screening key regulatory factors and signaling nodes from GRN are as follows: Use the filter function to screen ligand-receptor pairs that satisfy p<0.05 in the communication matrix constructed by CellChat, and extract the communication path name, probability, and significance. The frequency and average communication strength of each ligand-receptor pair were calculated using the summarise function, and the two were combined into a candidate gene pair score matrix. The screening threshold was set as the frequency being above the top 30 percentile and the average communication strength being above 0.

01. The screening results are the hub gene list. The key transcription factors obtained from RcisTarget analysis are then combined to identify key regulatory factors and signaling nodes. And / or, input the key regulatory factors and signal nodes into the drug-target database for retrieval, export the retrieval results, and visualize the undirected drug-target network using the igraph package; based on the drugs corresponding to the known targets included in the database, if the current indication of the drug is different from the disease being studied, then the repurposing of this drug can be considered in the future; if no known drug has been found for this target, then new drug development for this target can be considered in the future.

8. A screening system for drug targets with spatiotemporal specificity, characterized in that, The screening system is used to implement the screening method according to any one of claims 1-7, the screening system comprising: (1) Quantitative reconstruction of the spatial localization and functional pattern of cells in tissues, including: The data preprocessing module is used to preprocess the single-cell transcriptome sequencing data. The data preprocessing includes: quality control and normalization of the single-cell transcriptome sequencing data, data integration, principal component analysis dimensionality reduction and Louvain unsupervised clustering, and cell type annotation. The spatial positioning reconstruction module is used to reconstruct the spatial positioning of single cells. The spatial positioning reconstruction includes: constructing a spatial positioning gene set, non-negative matrix factorization, and reconstructing the spatial positioning of each cell subpopulation. The functional phenotype reconstruction module is used to reconstruct the functional patterns of single-cell biology. The functional phenotype reconstruction includes: single-sample gene set enrichment analysis, UMAP nonlinear dimensionality reduction, K nearest neighbor map construction and Leiden unsupervised clustering, and output of functional feature classification results and enrichment scores. (2) A spatiotemporally specific drug target discovery module, which includes: Cell-to-cell communication analysis module, used for performing cell-to-cell communication analysis; A module for constructing gene regulatory networks, used to build GRN gene regulatory networks centered on specific tissue microenvironment states; and, The target discovery module is used to screen key regulatory factors and signaling nodes from the gene regulatory network GRN. The obtained hub genes and their connected core nodes are potential spatiotemporally specific drug targets.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the screening method according to any one of claims 1-7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the screening method according to any one of claims 1-7.

Images