A Virtual Perturbation Analysis System for Targets Based on Single-Cell Transcriptome Data

Abstract

This application relates to the fields of bioinformatics and computer data processing technology, and discloses a target virtual perturbation analysis system based on single-cell transcriptome data, including a data input module, a spatial construction module, a virtual cell perturbation layer, a perturbation evaluation layer, and a result output module. The system constructs a cell state space by reducing the dimensionality of the expression matrix, and calculates the average direction vector for low-expression K-nearest neighbors of cells in the space to obtain the virtual perturbation displacement vector and virtual cell state. The system uses Wasserstein-1 distance to calculate the displacement of the functional module score distribution before and after perturbation to generate a direct perturbation score, and constructs a graphical Laplacian matrix based on the ligand-receptor expression parameters of the sending and receiving cells. The direct perturbation score is propagated in the intercellular communication network using a thermonuclear diffusion matrix, and aggregated to obtain a tissue-level reversal score. This system realizes the quantitative evaluation of the perturbation effects of target genes and candidate drugs from single-cell resolution to the tissue level.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and computer data processing technology, specifically to a target virtual perturbation analysis system based on single-cell transcriptome data. Background Technology

[0002] With the development of high-throughput sequencing technology, single-cell transcriptome sequencing has been widely used in the study of gene expression patterns and disease mechanisms. In drug development and disease target discovery, simulating the perturbation effect of reduced gene expression through computational methods is an important technique for evaluating candidate target genes and the results of drug intervention.

[0003] However, existing methods for calculating virtual gene perturbations mainly focus on analyzing expression changes at the single-cell level or within specific cell types, lacking a quantitative evaluation mechanism for networked propagation at the tissue level. Specifically, at single-cell resolution, existing methods struggle to generate perturbed virtual cell states through deterministic spatial coordinate translation calculations, leading to uncertainty in quantifying the single-cell perturbation intervention process. Furthermore, current technologies typically compare average gene expression values ​​before and after perturbation directly, failing to measure changes in functional module scores from a statistical distribution shift perspective, and thus unable to convert discrete single-cell state changes into objective evaluation values ​​for the entire cell type community. More importantly, biological tissues are holistic entities containing complex intercellular communication networks. Existing single-cell evaluation schemes assess intervention results for individual cell types in isolation, failing to construct a communication network matrix using ligand and receptor expression parameters between sending and receiving cells, and lacking a deterministic mathematical diffusion model to quantitatively propagate local direct perturbation scores along the intercellular communication network. Consequently, they cannot output quantitative evaluation values ​​that reflect the comprehensive effects of target genes or candidate drugs at the tissue level. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a target virtual perturbation analysis system based on single-cell transcriptome data, which solves the problem that existing technologies at single-cell resolution are unable to quantify and evaluate the comprehensive effects of perturbation at the tissue level, which is composed of intercellular communication networks.

[0005] To achieve the above objectives, the present invention provides the following technical solution: This invention provides a target virtual perturbation analysis system based on single-cell transcriptome data. The system receives a single-cell transcriptome expression matrix and a set of target genes, and performs virtual perturbation analysis and evaluation through the following data processing and algorithm execution flow: In the data preprocessing and spatial construction stage, the system performs quality control on the received single-cell transcriptome expression matrix, retaining cells with the number of detected genes within a preset range and the proportion of mitochondrial-derived reads below a preset threshold, while removing underexpressed genes. Subsequently, the system normalizes and logarithmically transforms the total expression of each cell in the single-cell transcriptome expression matrix, filters out a preset number of hypervariable genes, and performs scaling, centering, and variance standardization on their expression values. The system further performs principal component analysis for dimensionality reduction, extracting the first few principal components to construct the cell state space.

[0006] In the virtual cell perturbation calculation stage, the system calculates a virtual perturbation displacement vector for the target gene in the cell state space. Specifically, the system constructs a K-nearest neighbor set with a preset number of nearest neighbors for any current cell, and selects samples with target gene expression values ​​lower than that of the current cell to form a low-expression nearest neighbor set. The system obtains a direction vector by subtracting the coordinates of the current cell from the coordinates of each nearest neighbor cell in the low-expression nearest neighbor set, performs an equal-weighted average on this vector, and multiplies it by a perturbation intensity scaling factor to calculate the virtual perturbation displacement vector for a single gene. For combined perturbations containing two or more target genes, the system calculates the average of the virtual perturbation displacement vectors corresponding to each component gene to obtain the combined perturbation displacement vector. The original coordinates of the current cell are added to the corresponding displacement vector to generate the perturbated virtual cell state.

[0007] For invalid cells at the boundary in the cell state space, the system sets the following processing rules: when the low-expression nearest neighbor set of the current cell is empty, its corresponding virtual perturbation displacement vector is occupied by a zero vector, and the cell is excluded from participating in the construction of the subsequent pre-perturbation and post-perturbation module score distribution; if the number of valid perturbation cells in a certain cell type is lower than the preset minimum number, the system marks the target gene as a missing value or fills it with zero according to the preset rules as the reverse score of the perturbation evaluation score under the cell type.

[0008] In the perturbation evaluation and distribution displacement calculation stage, the system divides the virtual cell state into communities based on cell type annotation information. The system constructs the pre-perturbation functional module score distribution and the corresponding post-perturbation functional module score distribution, where the original cell functional module scores are pre-calculated using GSVA, AUCell algorithm, UCell algorithm, or AddModuleScore method. The system uses Wasserstein-1 distance to calculate the distribution displacement between the pre-perturbation and post-perturbation score distributions, and assigns positive or negative signs based on whether the direction of distribution change conforms to the preset treatment direction, generating a signed inverse score as the direct perturbation score at the cell type level.

[0009] In the tissue-level communication and propagation phase, the system extracts ligand and receptor gene expression parameters for each cell type based on the ligand-receptor relationship set. The system multiplies the average ligand expression value and ligand positivity rate of the sending cell type with the average receptor expression value and receptor positivity rate of the receiving cell type to calculate the communication score for the ligand-receptor pair. This is then summarized to obtain a directed communication adjacency matrix, which is symmetricized to construct an undirected weighted propagation matrix. Subsequently, the system calculates the degree matrix based on the undirected weighted propagation matrix, and then constructs a graph Laplacian matrix. The system uses this graph Laplacian matrix and preset thermonuclear diffusion parameters to construct a thermonuclear diffusion matrix, which is then multiplied by the cell type-level direct perturbation scoring matrix to calculate the post-propagation scoring matrix. Based on this, a tissue-level reversal score is obtained.

[0010] During the candidate drug results output phase, the system compares the target genes of the candidate drugs with the target gene set, and calls upon the tissue-level reversal scores corresponding to the single gene or combination of perturbation objects. The system converts the scores that match the treatment direction into positive therapeutic reversal intensity, sorts the candidate drugs from high to low according to the intensity value, and outputs the results; for candidate drugs that do not match the target genes, the system marks them as unscored and excludes them from the sorting or marks them as missing values.

[0011] This invention provides a target virtual perturbation analysis system based on single-cell transcriptome data. It has the following beneficial effects: 1. The invention achieves deterministic quantization of virtual cell state coordinates after target gene perturbation at single-cell resolution by screening low-expression nearest neighbor sets in the cell state space after principal component dimensionality reduction and calculating the average direction vector to obtain virtual perturbation displacement vector.

[0012] 2. This invention calculates the distribution displacement between the functional module score distribution before and after the perturbation based on Wasserstein-1 distance, and generates a signed inverse score accordingly, thereby transforming the discrete single-cell spatial coordinate translation into an objective evaluation value of the overall change of a specific cell type community.

[0013] 3. This invention utilizes the expression parameters of ligands and receptors in sending and receiving cells to construct a graph Laplacian matrix. By multiplying a preset thermonuclear diffusion matrix with a cell type-level direct perturbation score matrix, the direct perturbation score of local cell types can be quantitatively diffused along the intercellular communication network, ultimately outputting a quantified tissue-level reversal score. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the overall architecture of the ScalePert virtual screening system provided in an embodiment of the present invention; Figure 2This is a heatmap of virtual cell k-sensitivity analysis provided in an embodiment of the present invention; Figure 3 A heatmap of virtual cell sensitivity analysis provided in an embodiment of the present invention; Figure 4 This is a virtual cell single-gene perturbation heatmap provided in an embodiment of the present invention; Figure 5 This is a virtual cell combination perturbation heatmap provided in an embodiment of the present invention; Figure 6 A heatmap of virtual organoid β-sensitivity analysis provided in an embodiment of the present invention; Figure 7 Intercellular communication diagram provided for embodiments of the present invention; Figure 8 This is a heatmap of single-gene perturbation of virtual organoids provided in an embodiment of the present invention; Figure 9 This is a perturbation heatmap of virtual organoid combinations provided in an embodiment of the present invention; Figure 10 A tissue-level drug ranking heatmap provided for embodiments of the present invention. Detailed Implementation

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

[0016] See attached document Figure 1 An embodiment of the present invention provides a target virtual perturbation analysis system based on single-cell transcriptome data. This system uses single-cell RNA sequencing data or single-nuclear RNA sequencing data as its foundation. Through single-cell data preprocessing, cell state space construction, target gene virtual perturbation, functional module reversal scoring, inter-cell communication propagation, and drug target mapping, a computational process is formed from single-cell transcriptome data to the ranking of candidate therapeutic targets and candidate drugs.

[0017] This invention relates to a target virtual perturbation analysis system based on single-cell transcriptome data, which specifically includes: a data input module, a spatial construction module, a virtual cell perturbation layer, a perturbation evaluation layer, and a result output module.

[0018] The data input module is used to obtain single-cell transcriptome expression matrix and target gene set, as well as cell type annotation information, functional module gene set, ligand-receptor relationship set and candidate drug target information.

[0019] The space construction module is used to construct the cell state space based on the single-cell transcriptome expression matrix.

[0020] The virtual cell perturbation layer is used to perform virtual perturbation analysis on the target gene in the cell state space. It calculates the virtual perturbation displacement vector based on the neighborhood state of each cell in the cell state space to obtain the perturbed virtual cell state.

[0021] The perturbation evaluation layer is used to calculate the perturbation evaluation score of the target gene based on the changes in the virtual cell state compared to before the perturbation. To further evaluate the propagation of the perturbation effect in the tissue environment, the perturbation evaluation layer can be specifically represented as a virtual organoid perturbation layer. The virtual organoid perturbation layer is used to propagate the perturbation score based on the intercellular communication network and form a comprehensive evaluation result at the tissue level or organoid level.

[0022] The results output module is used to output evaluation results for the target gene based on the perturbation evaluation score.

[0023] Specifically, the data input module receives a single-cell or single-nuclear RNA sequencing expression matrix, and receives corresponding cell type annotation information, target gene set, functional module gene set, ligand-receptor relationship set, and candidate drug target information. The target gene set is used to define the gene objects to be subjected to virtual perturbation analysis; the functional module gene set is used to evaluate the impact of perturbation on disease-related functional states; the ligand-receptor relationship set is used to construct a communication network between cell types; and the candidate drug target information is used to map candidate drugs to single-gene perturbation objects or combined gene perturbation objects.

[0024] During system processing, the expression matrix is ​​transformed into a cell state space after quality control, standardization, and dimensionality reduction. Within this cell state space, the system determines the direction of virtual perturbation based on the target gene expression level and obtains the perturbed virtual cell state. Subsequently, based on cell type annotations, the system divides the virtual perturbation-post-perturbation states at the single-cell level into communities according to cell type, constructs the pre- and post-perturbation module score distributions for each cell type, calculates the reversal score, and obtains the direct perturbation score at the cell type level.

[0025] To further evaluate the propagation of perturbation effects within the tissue environment, the system constructs an intercellular communication network based on ligand-receptor relationships. Direct perturbation scores at the cell type level are then input into this network for propagation, yielding a tissue-level perturbation score. This tissue-level perturbation score is used to rank individual target genes, combinations of target genes, and candidate drugs.

[0026] Based on the above system architecture, the overall workflow of this invention includes the following steps.

[0027] S100 acquires single-cell or single-nuclear RNA sequencing expression matrices, and obtains cell type annotation information, target gene sets, functional module gene sets, ligand-receptor relationship sets, and candidate drug target information.

[0028] S200 preprocesses the expression matrix and constructs the cell state space based on the preprocessed expression matrix.

[0029] S300 performs virtual perturbation analysis on the target gene in the cell state space to obtain the perturbed virtual cell state.

[0030] S400 calculates the reversal score corresponding to the target gene based on the functional module score distribution before and after the perturbation.

[0031] S500: When the target perturbation object includes two or more target genes, the perturbation displacement vectors corresponding to each target gene are combined in the cell state space, and the combined perturbation score is calculated based on the combined displacement.

[0032] S600, based on cell type annotation, divides the virtual perturbation state at the single-cell level into communities according to cell type, constructs the pre- and post-perturbation module score distribution at the cell type level, and thus calculates the direct perturbation score at the cell type level.

[0033] S700 constructs an intercellular communication network based on ligand-receptor relationships and propagates direct perturbation scores on the intercellular communication network to obtain tissue-level perturbation scores.

[0034] S800 maps candidate drug targets to single-gene perturbation objects or combined gene perturbation objects, and outputs candidate drug rankings based on tissue-level perturbation scores.

[0035] S900 outputs the ranking of candidate therapeutic targets, the ranking of candidate target combinations, the ranking of candidate drugs, and the corresponding visualization results.

[0036] In the above process, data input and preprocessing, cell state space construction, virtual perturbation analysis, reversal score calculation, communication network propagation, and drug target mapping are sequentially linked. The output of the previous processing step serves as the input for the next processing step, enabling the system to obtain virtual perturbation evaluation results at the cell type and tissue levels from single-cell transcriptome data.

[0037] Therefore, this invention establishes a target virtual perturbation analysis system based on single-cell transcriptome data. This system achieves tissue-level evaluation of target genes, target gene combinations, and candidate drugs through virtual gene perturbation calculations in the cell state space and perturbation propagation calculations in the intercellular communication network.

[0038] Within the overall system architecture described above, this invention first performs unified processing on single-cell transcriptome data to obtain the cell state space required for subsequent virtual perturbation calculations. Single-cell transcriptome data includes single-cell RNA sequencing data and single-nuclear RNA sequencing data, both using the cell or cell nucleus as the basic unit of analysis. For ease of description, they are collectively referred to as single-cell transcriptome data below.

[0039] The system receives a single-cell transcriptome expression matrix as its basic input. The expression matrix is ​​arranged with cells as rows and genes as columns, and each matrix element represents the expression value of the corresponding gene in that cell. This expression matrix is ​​denoted as: in, This represents the input single-cell transcriptome expression matrix. Indicates the number of cells. Indicates the number of genes. Indicates the first The first cell The expression values ​​of each gene.

[0040] The expression values ​​include raw count values ​​or expression values ​​after basic processing. The system converts them into expression data suitable for dimensionality reduction and nearest neighbor analysis based on a unified data processing flow.

[0041] In addition to the expression matrix, the system receives cell type annotation information corresponding to the expression matrix. This cell type annotation information identifies the cell population to which each cell belongs and corresponds to the cell number in the expression matrix. Subsequently, the system uses this annotation information to classify the perturbation states at the single-cell level to the corresponding cell populations, in order to construct a module score distribution at the cell type level.

[0042] The system also receives a target gene set, a functional module gene set, a ligand-receptor relationship set, and candidate drug target information. The target gene set is used to identify the gene objects that need to be subjected to virtual perturbation analysis; the functional module gene set is used to calculate the functional module score of the cell; the ligand-receptor relationship set is used to construct the communication network between cell types; and the candidate drug target information is used to map candidate drugs to single-gene perturbation objects or combined gene perturbation objects.

[0043] The above input data is matched with the expression matrix at the gene name or gene number level. Genes that cannot be matched with the expression matrix are not included in the corresponding calculation.

[0044] To reduce the impact of low-quality data on the cell state space, the system performs quality control on the cells and genes in the input expression matrix. In one specific implementation, the system retains cells with the number of detected genes within a preset range and the proportion of mitochondrial-derived reads below a preset threshold, and removes genes that are underexpressed in the cell population.

[0045] The specific thresholds for the number of genes to be detected, the proportion of mitochondria, and the filtering of low-expression genes can be set by those skilled in the art based on the quality of single-cell transcriptome data and the sample type.

[0046] The expression matrix after quality control undergoes standardization. The system normalizes the total expression level for each cell to ensure comparable sequencing depths across different cells. Subsequently, a logarithmic transformation is applied to the normalized expression values ​​to reduce the impact of highly expressed genes on subsequent dimensionality reduction results. The expression matrix after total expression normalization and logarithmic transformation is then used for screening for highly variable genes.

[0047] The system screens for highly variable genes based on a standardized expression matrix. Highly variable genes are those exhibiting significant expression variation across cell populations, used to retain expression information relevant to differences in cell state. In one specific implementation, the system selects a predetermined number of highly variable genes as input for PCA dimensionality reduction. After screening, the system scales the expression values ​​of the selected genes to ensure a consistent numerical scale across different genes in the dimensionality reduction analysis.

[0048] For expression value scaling, centering, and variance normalization, those skilled in the art can use conventional methods in the field of single-cell transcriptome analysis.

[0049] When input data comes from multiple samples or multiple sequencing batches, the system corrects for technical differences between samples. Batch correction is used to reduce the impact of differences in sample source, sequencing batch, or platform on the construction of the cell state space, so that the subsequent PCA space mainly reflects the differences in cell transcriptional state. Batch correction can be performed at the expression matrix level or at the dimensionality reduction representation level, and the specific tools and parameters do not limit the scope of protection of this invention.

[0050] After the above preprocessing, the system performs principal component analysis based on the processed expression matrix to obtain the cell embedding matrix in the PCA space: ; in, This represents the embedding matrix of cells in the PCA space. Indicates the number of cells. This indicates the spatial dimension of PCA.

[0051] No. The coordinates of each cell in PCA space are: ; in, Indicates the first Low-dimensional transcriptome state representation of individual cells.

[0052] This PCA space does not directly represent the expression value of a single gene, but rather the cellular transcriptome state formed by multiple genes. The system performs nearest neighbor search in this space, enabling subsequent virtual perturbations to be calculated based on local cellular state relationships.

[0053] In one specific implementation, the system performs total normalization and logarithmic transformation on single-cell transcriptome data, screens for highly variable genes, scales expression values, and extracts the first few principal components as the cell state space. This cell state space serves as the common basis for subsequent target gene expression vector extraction, K-nearest neighbor construction, low-expression neighbor screening, and perturbation shift calculation.

[0054] After the above processing, the system obtains standardized expression data, cell type annotations, target gene expression data, functional module gene sets, and PCA embedding matrices. These data are then fed into the subsequent virtual cell perturbation process to calculate the perturbation displacement in the target gene depletion direction within the same cell state space.

[0055] See attached document Figure 4 The system performs virtual perturbation analysis on the target gene in this PCA space and calculates a reversal score based on the changes in the functional module score distribution before and after the perturbation. This process is used to evaluate the impact of reduced target gene expression on disease-related functional status.

[0056] In this invention, the virtual perturbation of the target gene does not directly modify the gene expression values ​​in the original expression matrix. Instead, it calculates the direction of cell movement based on the local cell state relationships in the PCA space. Each point in the PCA space corresponds to the transcriptomic state of a cell, and cells with low target gene expression located in the current cell's neighborhood provide a local reference after the target gene expression has decreased. In this way, the system can construct a virtual cell state after target gene depletion in observational single-cell transcriptomic data.

[0057] For the target gene to be perturbed The system extracts the expression vector of the gene in all cells from the standardized expression data: ; in, Indicates target gene Expression vector in all cells, Indicates the number of cells. Indicates target gene In the Expression values ​​in individual cells.

[0058] This expression vector corresponds one-to-one with the cell number in the aforementioned PCA embedding matrix.

[0059] The system constructs a K-nearest neighbor set for the i-th cell in the PCA space, denoted as . .

[0060] in, Represents the nearest neighbor number. Indicates the relationship with the first The set of the K nearest neighbors of a given cell in the PCA space. For the distance metric of nearest neighbor search in the PCA space, those skilled in the art can use conventional distance metrics such as Euclidean distance and cosine distance.

[0061] After obtaining the K-nearest neighbor set, the system determines the target gene. The expression value is used to filter the low-expression nearest neighbors in the K-nearest neighbor set: ; in, Indicates the first Each cell corresponds to the target gene. Low-expression nearest neighbor set; Indicates the cell number of the nearest neighbor; Indicates target gene In neighboring cells The expression value in; Indicates target gene In the The expression values ​​in each cell. Cells in this set simultaneously satisfy two conditions: Firstly, it belongs to the first The nearest neighbors of each cell in PCA space; Secondly, target genes The expression value is lower than the first Each cell.

[0062] After determining the low-expression nearest neighbor set, the system uses the average direction of cells in this set relative to the current cell as the local movement direction after the target gene expression is reduced. For cells whose low-expression nearest neighbor set is not empty, the target gene... For the first The virtual perturbation displacement vector generated by each cell is: ; in, Indicates target gene For the first The virtual perturbation displacement vector generated by each cell; This represents the disturbance intensity scaling factor; This indicates the number of cells in the low-expression nearest neighbor set; Indicates the first The coordinates of each cell in PCA space; Indicates low expression in neighboring cells Coordinates in PCA space; Indicates by the first The direction vector of a cell pointing to its low-expression neighboring cell.

[0063] The above formula shows that the first The direction of perturbation in each cell is determined by the cells in its local neighborhood that express the target gene at low levels. (Scaling factor) Used to control the amplitude of the disturbance displacement, enabling the system to adjust the virtual exhaustion level in the same local direction.

[0064] For cells with an empty nearest neighbor set for low expression, the system marks these cells as invalid perturbation cells for the target gene and places their perturbation shift as a zero vector. Invalid perturbation cells do not participate in the construction of the pre-perturbation module score distribution and post-perturbation module score distribution corresponding to the target gene. If the number of valid perturbation cells in a certain cell type is lower than a preset minimum number, the system marks the reverse score of the target gene in that cell type as a missing value, or fills it with zeros according to preset rules.

[0065] Based on the perturbation displacement vector, the system calculates the target gene. After the disturbance The virtual position of each cell in PCA space: ; in, Indicates target gene After virtual perturbation The system calculates the position of each cell in the PCA space. It then performs the above calculations on each cell to obtain the target gene. The set of virtual cell states after perturbation.

[0066] Since the perturbation location lies within the PCA space, the system needs to convert this location into a functional module score to evaluate the impact of the perturbation on disease-related functional states. The system pre-calculates the module score for each primitive cell based on the functional module gene set. A functional module is a set of genes consisting of multiple function-related genes used to characterize disease-related biological processes.

[0067] In one specific implementation, the functional modules include an innate immune sensing module, an oxidative metabolic stress module, a DNA damage repair module, and a DNA sensing module.

[0068] The module score for each primitive cell can be calculated using conventional algorithms in the field of single-cell transcriptomics, such as, but not limited to, GSVA (Gene Set Variation Analysis), AUCell algorithm, UCell algorithm, or the AddModuleScore method of the Seurat platform. For the perturbation position, the system searches for a preset number of original cells that are closest to the perturbation position in the original PCA space, and reads the functional module scores pre-calculated by the original cells; In one specific implementation, the system uses the average score of the functional modules of the nearest primitive cells as the module score of the perturbed location, or uses a weighted average that is inversely proportional to the distance as the module score of the perturbed location. Thus, the perturbed PCA spatial location can be mapped to the corresponding functional module score.

[0069] For nearest neighbor interpolation estimation, those skilled in the art can use the K-nearest neighbor interpolation method.

[0070] For functional modules The system constructs a pre-perturbation module score distribution based on the module scores of effective cells before the perturbation, denoted as [the distribution is missing from the original text]. ; The module scores corresponding to the perturbed virtual cell states form the perturbed module score distribution, denoted as [distribution]. .in, Indicates the first Each functional module This represents the distribution of module scores before the disturbance. This represents the distribution of module scores after the perturbation.

[0071] When cell type level scores are required, the pre-perturbation module score distribution and the post-perturbation module score distribution are constructed separately according to cell type. When a global score needs to be output, the pre-perturbation module score distribution and the post-perturbation module score distribution can be constructed based on all effective perturbation cells.

[0072] Preferably, the direct perturbation score entering the intercellular communication propagation process is a cell type-level direct perturbation score calculated separately for each cell type.

[0073] The system uses the Wasserstein-1 distance to calculate the difference between the module score distributions before and after the disturbance, obtaining a non-negative distributed displacement. The Wasserstein-1 distance is used to measure the overall displacement between the score distributions of two modules; therefore, the system uses the distribution change as the evaluation object of the disturbance effect, rather than just a single statistic.

[0074] After obtaining the non-negative distributed displacement, the system generates a signed inverse score based on the direction of change of the module score distribution after perturbation relative to the module score distribution before perturbation. When the preset treatment direction is a decrease in module score, and the module score distribution after perturbation is lower than the module score distribution before perturbation, the system assigns the distributed displacement a treatment direction sign; when the module score distribution after perturbation increases or does not conform to the preset treatment direction, the system assigns the distributed displacement a non-treatment direction sign. Thus, the signed inverse score simultaneously contains information on the magnitude and direction of the distributed displacement before and after perturbation.

[0075] When ranking candidate therapeutic targets, combinations of candidate targets, and candidate drugs, the system converts signed reversal scores that align with the therapeutic direction into positive therapeutic reversal strength and ranks them based on the therapeutic reversal strength.

[0076] In one specific implementation, the system screens effective perturbation cells before calculating the reversal score. The system retains cells whose perturbation displacement amplitude meets preset conditions and constructs a pre- and post-perturbation module score distribution based on the retained cells. This process reduces the impact of cells without effective displacement on the score calculation.

[0077] To evaluate the statistical support of the target gene perturbation score, the system permutes the correspondence between module scores and cells, constructing a random null distribution. In each permutation, the system keeps the cell's position in the PCA space and the target gene expression vector unchanged, shuffles the correspondence between module scores and cells, and recalculates the perturbation score.

[0078] The system obtains the statistical significance of the perturbation score by comparing it with a random null distribution. For test results generated from simultaneous multi-gene, multi-module, and multi-cell type analyses, the system performs Benjamini-Hochberg FDR correction.

[0079] See attached document Figure 2 and attached Figure 3 In the virtual cell perturbation analysis, the system varied the nearest neighbor number and the perturbation intensity scaling factor, and repeatedly performed virtual perturbation and reversal score calculations for the target gene under different parameter settings. The system compared the consistency of the target gene rankings obtained under different parameter settings to evaluate the dependence of the virtual cell perturbation results on the parameter settings.

[0080] Through the above processing, the system can output a virtual perturbation reversal score for each target gene in different functional modules and cell types. This result serves as input for subsequent combined perturbation analysis and tissue-level communication propagation analysis.

[0081] See attached document Figure 5After completing the calculation of single-gene virtual perturbation and reversal scores, the system further analyzes the combined perturbation of two or more target genes. The combined perturbation is used to characterize the impact of the simultaneous reduction of expression of multiple target genes on cell state and functional module scores, and its calculation basis is still the aforementioned single-gene perturbation shift vector.

[0082] In this invention, the combined perturbation does not reconstruct independent expression matrices, but rather reuses the perturbation shifts of each target gene within the same PCA cell state space. The system first obtains the single-gene perturbation shifts corresponding to each component gene, and then superimposes these perturbation shifts in an equal-weighted manner. Thus, the combined perturbation and the single-gene perturbation maintain the same cell state space and scoring system, allowing for direct entry into the module score estimation and reverse score calculation processes.

[0083] In combined perturbation, if a component gene does not have effective low-expression neighbors in a certain cell, the system can record the perturbation displacement vector of that component gene in that cell as a zero vector for combined calculation, or calculate the combined perturbation displacement only based on the component genes with effective perturbation displacement in that cell. Preferably, the system uses a zero vector placeholder and, after the combined perturbation displacement calculation is completed, screens effective combined perturbation cells based on the magnitude of the combined perturbation displacement.

[0084] Let the target gene set corresponding to the combined perturbation be . ,in It contains two or more target genes. For the first The system has obtained a set of cells based on the aforementioned virtual cell perturbation process. Each target gene The corresponding disturbance displacement vector When a certain target gene In the When there are no effectively low-expression neighbors in a cell, the system will perturb the shift vector corresponding to the target gene. Denote it as the zero vector, or denote the target gene from the first... Remove from the effective combination of components of each cell.

[0085] In one specific implementation, the system employs a zero-vector placeholder approach and defines the combined perturbation displacement vector as the average value of the perturbation displacement vectors of each component gene: ; in, Represents the target gene set; Represents the target gene set The number of target genes; This represents the dimension of the PCA cell state space.

[0086] When target gene In the When there are no effectively low-expression nearest neighbors in a cell, the system records the single-gene perturbation shift vector corresponding to the target gene as the zero vector: ; After the combined perturbation, the first The virtual positions of each cell in PCA space are: ; in, Represents the target gene set After combined perturbation The position of each cell in PCA space; Indicates the first The coordinates of each cell in the PCA space before perturbation.

[0087] After obtaining the virtual locations after combined perturbations, the system estimates the perturbed functional module scores in the same manner as with single-gene perturbations. Specifically, the system searches for and matches the original PCA space. The module score after combined perturbation is estimated based on the functional module scores of neighboring primitive cells. For K-nearest neighbor search and score interpolation, those skilled in the art can use conventional implementation methods in single-cell data analysis.

[0088] For functional modules The system constructs a pre-perturbation score distribution from the module scores before the perturbation and a post-perturbation score distribution from the module scores after the perturbation. The Wasserstein-1 distance is then used to calculate the distribution shift between the two distributions. This yields the target gene set. For functional modules The combined perturbation reversal score.

[0089] This score has the same meaning as the single-gene perturbation reversal score, both used to characterize the changes in the distribution of functional module scores before and after the perturbation.

[0090] When calculating the combined perturbation reversal score, the system selects effective combined perturbation cells based on the combined perturbation displacement amplitude. For cells whose combined perturbation displacement amplitude does not meet the preset conditions, the system does not include them in the construction of the pre-combined perturbation score distribution and the post-combined perturbation score distribution. If the number of effective combined perturbation cells in a certain cell type is lower than the preset minimum number, the system marks the combined perturbation reversal score of the target gene set in that cell type as a missing value, or fills it with zero according to preset rules.

[0091] To enable combined perturbation results at the single-cell level to enter the subsequent tissue-level propagation process, the system aggregates perturbation scores at the cell type level based on cell type annotations. Cell type annotations correspond to cell numbers in the expression matrix, thus allowing the system to determine the cell type of each valid perturbation cell and aggregate perturbation results within the same cell type into a direct perturbation score for that cell type.

[0092] In one specific implementation, the system establishes scoring results using the target perturbation object, functional module, and cell type as indexes. The target perturbation object includes a single target gene and a set of target genes; the functional module is a preset disease-related functional module; and the cell type is the annotated cell population in the input single-cell transcriptome data.

[0093] To obtain direct perturbation scores at the cell type level, the system groups cells by target perturbation object, functional module, and cell type. For any target perturbation object, any functional module, and any cell type, the system selects cells belonging to that cell type that meet the valid perturbation conditions, constructs the pre-perturbation module score distribution and the post-perturbation module score distribution within that cell type, and calculates the signed inverse score of the target perturbation object under that cell type and that functional module based on the two distributions.

[0094] Subsequently, the system summarizes the signed inversion scores of multiple functional modules according to the preset module combination rules to obtain the cell type level direct perturbation score of the target perturbation object under the cell type.

[0095] Cell type-level direct perturbation score refers to the direct impact of a target perturbation object on the functional module score distribution within a specific cell type. This score has not yet been propagated through the inter-cell communication network, and therefore is used to characterize the direct computational result of the target perturbation object in the corresponding cell type. The subsequent thermonuclear diffusion propagation module uses this score as input to propagate it along the inter-cell communication network to the tissue level.

[0096] The system organizes the direct perturbation scores of multiple cell types into a scoring matrix, denoted as . ; in, This represents the cell type level direct perturbation score matrix. Rows in this matrix correspond to cell types, and columns correspond to the target perturbation objects.

[0097] Before forming the scoring matrix, the system combines the inverse scores of multiple functional modules of the same target perturbation object under the same cell type into a comprehensive direct perturbation score according to preset module combination rules. Module combination rules include equal-weighted averaging, preset-weighted averaging, averaging only the scores of the treatment direction modules, or weighted aggregation based on disease relevance weights. The target perturbation objects include single-gene perturbation objects, combined-gene perturbation objects, and perturbation objects mapped from candidate drug targets. The elements in the matrix represent the direct perturbation scores of the corresponding target perturbation object in the corresponding cell type.

[0098] During cell type-level aggregation, the system preserves the directional attributes of the perturbation scores. Treatment direction reversal and non-treatment direction changes are recorded in the cell type-level direct perturbation score matrix according to a unified sign rule, enabling subsequent tissue-level propagation to be calculated using the same score matrix. For ranking, the system further converts the propagated signed scores into positive therapeutic reversal intensity. For scores of different functional modules, the system forms a comprehensive score for the target perturbation object according to preset module combination rules; these rules can employ averaging or weighted summarization, aiming to convert perturbation effects across multiple functional modules into a unified cell type-level input.

[0099] See attached document Figure 5 The combined perturbation results can be displayed in the form of a heatmap, where the horizontal axis corresponds to the target gene combination and the vertical axis corresponds to the cell type. The color or value represents the reversal score of the combined perturbation in the corresponding cell type. This result is used to show the perturbation effect of different target gene combinations in various cell types and serves as input for subsequent virtual organoid perturbation analysis.

[0100] Through the aforementioned combined perturbation and cell type-level aggregation processing, the system organizes the displacement calculations and module scoring results at the single-cell level into a cell type-level direct perturbation scoring matrix. This matrix preserves both the differences between different target perturbation objects and the response differences between different cell types, thus enabling it to connect with the inter-cell communication network construction module.

[0101] See attached document Figure 7 After obtaining the direct perturbation score at the cell type level, the system further constructs a communication network based on intercellular ligand-receptor relationships and propagates the perturbation score on this network. Through this process, the direct perturbation results at the single cell type level are converted into tissue-level perturbation results that include intercellular communication relationships.

[0102] The input to the intercellular communication network includes cell type annotation information, gene expression data for each cell type, and a predefined set of ligand-receptor relationships. The ligand-receptor relationship set includes ligand genes and their corresponding receptor genes, used to describe the signal transduction relationships between different cell types. In one specific embodiment, the ligand-receptor relationship set includes CCL5-CCR5, CCL2-CCR2, CXCL10-CXCR3, IL1B-IL1R1, IL18-IL18R1, TNF-TNFRSF1A, IFNG-IFNGR1, IL6-IL6R, HMGB1-TLR2, HMGB1-TLR4, HMGB1-RAGE, VCAM1-ITGA4, ICAM1-ITGAL, VEGFA-FLT1, VEGFA-KDR, IL1A-IL1R1, CXCL8-CXCR1, and CXCL8-CXCR2.

[0103] The system groups single-cell expression data by cell type based on cell type annotation information. For any sending and receiving cell type, the system reads the expression level of the ligand gene in the sending cell type and the expression level of the corresponding receptor gene in the receiving cell type. Based on ligand and receptor expression, the system calculates the communication strength between the sending and receiving cell types.

[0104] In one specific implementation, for the sending cell type and the receiving cell type, the system calculates the average normalized expression value of the ligand gene and the proportion of positive cells in the sending cell type for each preset ligand-receptor pair, and calculates the average normalized expression value of the receptor gene and the proportion of positive cells in the receiving cell type. The average ligand expression value, ligand positivity rate, average receptor expression value, and receptor positivity rate were then multiplied to obtain the ligand-receptor pair in... to The communication score in the direction. The system sums or weighted sums the communication scores of all ligand-receptor pairs to obtain the communication strength from the sending cell type to the receiving cell type.

[0105] To reduce the impact of random expression fluctuations on the communication network, the system performs a permutation test on the communication relationships between cell types. In one specific implementation, the system shuffles the cell type label or ligand-receptor expression correspondences and recalculates the communication strength between cell types, thereby forming a random background distribution. The system compares the true communication strength with the random background distribution and retains the communication relationships that meet the preset significance criteria. The filtered communication relationships are used to construct the inter-cell type communication adjacency matrix.

[0106] The adjacency matrix for inter-cell type communication is denoted as: ; in, Represents the adjacency matrix for inter-cell type communication. Indicates the number of cell types. Indicates the first The cell type and the first The communication strength between cell types. Rows in the matrix correspond to sending cell types, and columns correspond to receiving cell types.

[0107] In one specific implementation, the system uses a communication adjacency matrix. To ensure that subsequent thermonuclear diffusion has a clear propagation meaning in cell-type networks, the system performs row normalization processing on the directed communication adjacency matrix. Converted into a nonnegative propagation matrix for thermonuclear diffusion Propagation matrix It can be by It can be obtained through normalization, or it can be obtained from... It is obtained by symmetrizing its transpose matrix.

[0108] Preferably, the system adopts We obtain the undirected weighted propagation matrix, where, for The transpose of the matrix, the undirected weighted propagation matrix It inherently possesses numerical symmetry. The system directly uses the undirected weighted propagation matrix. Construct the degree matrix .

[0109] The system is based on the propagation matrix Construct the degree matrix Degree matrix Let be a diagonal matrix, and its first... The diagonal elements are the propagation matrix. The Middle The sum of connection weights corresponding to each cell type. Based on the propagation matrix. Sum-degree matrix The system constructs the graph Laplace matrix: ; The Graph Laplace matrix is ​​used to represent the topology of cell-type communication networks and serves as the fundamental matrix for thermonuclear diffusion propagation.

[0110] The cell type level direct perturbation score matrix is ​​denoted as ,in The rows correspond to cell types, and the columns correspond to target perturbation objects. The system calculates the thermonuclear diffusion matrix. and through Obtain the perturbation score matrix after propagation, where, This is a thermonuclear diffusion parameter used to control the diffusion scale of perturbation scores in intercellular communication networks.

[0111] Therefore, thermonuclear diffusion propagates perturbation scores along the cell type dimension, resulting in a matrix after propagation. The cells are still grouped by cell type and the target perturbation object by column.

[0112] The aforementioned thermonuclear diffusion propagation is used to propagate the direct perturbation score along the cell type communication network. If a communication connection exists between two cell types, the direct perturbation score propagates between the relevant cell types according to the communication network structure; if no effective communication connection is formed between two cell types, the propagation relationship is defined by the network structure in the Graph Laplace matrix. Thus, the propagated score... It includes not only direct perturbation information in the original cell type, but also propagation information introduced by intercellular communication relationships.

[0113] See attached document Figure 8 For single-gene perturbation targets, the system inputs their cell type-level direct perturbation scores into the aforementioned thermonuclear diffusion propagation process to obtain post-propagation scores for each cell type, which are then further aggregated into tissue-level reversal scores. These tissue-level reversal scores are used to rank different single-gene perturbation targets.

[0114] See attached document Figure 9 For combined perturbation objects, the system also inputs the cell type level direct perturbation score of the combined perturbation into the intercellular communication network, and obtains the tissue level reverse score of the combined perturbation after thermonuclear diffusion propagation. Since both the combined perturbation score and the single gene perturbation score use the same cell type index and scoring direction rules, they can enter the same propagation process for calculation.

[0115] In one specific implementation, the system converts the propagated signed perturbation score into therapeutic reversal intensity and then performs tissue-level aggregation.

[0116] Specifically, for the same target perturbation object, the system reads its post-propagation score across different cell types, converts scores that align with the treatment direction into positive reversal intensity, and generates a tissue-level reversal score using averaging, summation, or a weighted average based on cell type proportions. A higher tissue-level reversal score indicates a stronger overall therapeutic reversal effect of the target perturbation object within the context of intercellular communication networks.

[0117] During aggregation, the system reads The propagation score of the same target perturbation object across different cell types is calculated and then aggregated according to a predefined aggregation rule to form a tissue-level score. This aggregation rule can be implemented by averaging, summing, or weighting by cell type. The resulting tissue-level inverse score is used to represent the overall perturbation effect of the target perturbation object in the context of a cell communication network.

[0118] See attached document Figure 6 The system repeatedly performs thermonuclear diffusion propagation and tissue-level scoring calculations under multiple preset thermonuclear diffusion parameters, and compares the consistency of target perturbation object ranking under different parameter settings. This process is used to evaluate the dependence of tissue-level propagation results on diffusion parameter settings.

[0119] Through the above processing, the system converts direct perturbation scores at the cell type level into inverse scores at the tissue level. This process allows the virtual perturbation results to go beyond direct scores within a single cell type, instead combining the communication networks between cell types to form tissue-level evaluation results, providing input for subsequent virtual screening of candidate drugs.

[0120] See attached document Figure 10 After obtaining tissue-level reversal scores for single-gene perturbations and combined perturbations, the system further receives candidate drug and target information, and converts the candidate drugs into target perturbation objects that can be calculated by the virtual perturbation analysis system. Through this process, the evaluation of candidate drugs is placed within the same computational framework as the aforementioned target gene perturbation analysis.

[0121] The system receives candidate drug target information including the candidate drug name and one or more target genes corresponding to that candidate drug. The target genes are matched against the target gene set. The matching results are used to determine the perturbation object corresponding to the candidate drug.

[0122] When there is one matching target gene in the target gene of a candidate drug, the system treats the candidate drug as a single gene perturbation object; when there are two or more matching target genes in the target gene of a candidate drug, the system treats the candidate drug as a combined perturbation object. When the target genes of a candidate drug fail to match any genes in the target gene set or expression matrix, the system marks the candidate drug as an unscored candidate drug and excludes it from the candidate drug ranking, or marks the candidate drug as a missing value for tissue-level reverse scoring.

[0123] In this invention, drug target mapping does not directly predict the binding relationship between drugs and molecules. Instead, it uses known drug target information to associate candidate drugs with the aforementioned virtual perturbation process. Therefore, drug evaluation can reuse cell state changes, functional module reversal scores, and tissue-level propagation results obtained from single-gene perturbations and combined perturbations.

[0124] Let the first The set of matchable target genes corresponding to each candidate drug is: .when When the candidate drug contains two or more target genes, the system calculates the perturbation shift corresponding to the candidate drug using a combined perturbation approach: ; in, Indicates the first The set of target genes corresponding to each candidate drug For the first The perturbation displacement vector generated by each cell; Represents the set of target genes The number of target genes that can be matched.

[0125] when When only one matching target gene is present, the system directly retrieves the single-gene perturbation shift and inversion score results corresponding to that target gene. Therefore, both single-target and multi-target drug candidates can enter the same scoring process.

[0126] when When it is an empty set, it represents the first set. If a candidate drug does not have a matching target gene that can enter the current virtual perturbation process, the system will mark the candidate drug as an unscored candidate drug and exclude it from the candidate drug ranking, or mark its cell type level direct perturbation score and tissue level inverse score as missing values.

[0127] The system calculates the perturbation-induced cell position based on the perturbation displacement corresponding to the candidate drug, and obtains the perturbation-induced module score distribution according to the aforementioned functional module score estimation method. Subsequently, the system divides the virtual cell state after perturbation of the candidate drug according to cell type annotation, constructs the functional module score distribution before and after perturbation within each cell type, and calculates the direct perturbation score of the candidate drug at each cell type level based on Wasserstein-1 distance.

[0128] For each candidate drug, the system inputs its cell type-level direct perturbation score into the intercellular communication network, and obtains the propagated perturbation score through thermonuclear diffusion. The propagated score is then aggregated according to a preset tissue-level aggregation rule to form the tissue-level reversal score of the candidate drug. The tissue-level reversal score is used to represent the comprehensive degree of reversal of disease-related functional modules by the candidate drug in the context of the intercellular communication network.

[0129] After calculating the tissue-level reversal scores for all candidate drugs, the system converts tissue-level reversal scores that align with the therapeutic direction into positive therapeutic reversal strengths and ranks the candidate drugs from highest to lowest based on these strengths. The ranking results are used to determine the priority of candidate drugs within the tissue context corresponding to the current single-cell transcriptome data. When outputting the ranking results, the system retains the drug name, corresponding target gene, cell type-level score, and tissue-level score, allowing users to simultaneously view the overall drug ranking and cell type origin.

[0130] In one specific implementation, the candidate drug target mapping includes 25 candidate drugs and their corresponding targets. The system maps each candidate drug to one or more hub gene perturbation objects and calculates its cell type-level direct perturbation score and tissue-level reversal score, respectively. (See Appendix) Figure 10 The system displays the scores of candidate drugs in different cell types and the ranking results at the tissue level using heatmaps and sorted bar charts.

[0131] The system outputs results including candidate therapeutic target ranking, candidate target combination ranking, and candidate drug ranking. Candidate therapeutic target ranking is derived from single-gene virtual perturbation results; candidate target combination ranking is derived from combined perturbation results; and candidate drug ranking is derived from tissue-level inverse scoring after drug target mapping. These three types of results use the same cell state space, functional module scoring system, and intercellular communication propagation process, thus enabling comparison and display within the same system.

[0132] To facilitate result verification, the system also outputs sensitivity analysis results for each parameter setting. (See attached document.) Figure 2 and attached Figure 3 The system outputs the sorting consistency results corresponding to the virtual cell perturbation parameters. (See attached document.) Figure 6 The system outputs the ranking consistency results corresponding to the thermonuclear diffusion parameters. Sensitivity analysis results are output along with the ranking results of candidate targets and candidate drugs to demonstrate the consistency of the system's calculations within the preset parameter range.

[0133] This invention also provides a computer implementation of the aforementioned virtual perturbation analysis system. The system is deployed in a computing device, which includes a processor and a memory. The memory stores a computer program, and when the processor executes the computer program, it performs single-cell transcriptome data input, PCA cell state space construction, virtual cell perturbation, Wasserstein reversal scoring, combined perturbation, intercellular communication network construction, thermonuclear diffusion propagation, and virtual drug screening.

[0134] The computer program can also be stored in a computer-readable storage medium. When the computer program is executed by a processor, it implements the virtual perturbation method based on single-cell transcriptome data of the present invention. The computer-readable storage medium includes a non-volatile storage medium capable of storing program code.

[0135] Through the aforementioned drug target mapping and result output process, the system incorporates candidate drug evaluation into a virtual perturbation calculation framework, enabling candidate drugs to be evaluated at the tissue level according to the perturbation effects of single genes or combinations of genes corresponding to their targets. Thus, the system can output complete virtual screening results from candidate targets, combinations of candidate targets to candidate drugs.

[0136] Specific application examples: The following is in conjunction with the appendix Figure 1 ~Attached Figure 10 The specific application process of the present invention will be described below. It should be understood that the following embodiments are only used to assist in illustrating the technical solution and application effects of the present invention, and are not intended to limit the scope of protection of the present invention.

[0137] Example 1: Overall Implementation of a Virtual Perturbation Analysis System Based on Single-Nuclear RNA Sequencing Data from Kidney Transplant Fibrosis This embodiment illustrates the overall implementation process of the virtual disturbance analysis system of the present invention, which can be referred to in conjunction with the appendix. Figure 1 Please provide an explanation.

[0138] In this embodiment, the input data was the single-nuclear RNA sequencing dataset GSE195718 from the NCBIGEO database. This dataset included 6 kidney transplant fibrosis samples and 3 normal kidney transplant biopsy samples. After Seuratv4 standard quality control and Harmony batch integration, high-quality single-cell / monocyte transcriptome data for 11 kidney cell types were obtained.

[0139] The system first preprocesses the input expression matrix. Specifically, this includes: normalizing the total expression values, performing a log1p transformation on the normalized expression values, screening the top 3000 hypervariable genes using the Seurat_v3 method, scaling the expression values ​​of the hypervariable genes, and extracting the top 50 principal components as the cell state space. This cell state space is denoted as... .

[0140] In this embodiment, the target gene set includes nine hub genes: TLR2, ITGAM, CCL5, HIF1A, VCAM1, JUN, AIM2, IFI16, and PYCARD. Among them, TLR2, ITGAM, and CCL5 are genes related to innate immune perception; HIF1A, VCAM1, and JUN are genes related to oxidative metabolic stress; and AIM2, IFI16, and PYCARD are genes related to DNA perception.

[0141] The system further defines four disease-related functional modules, namely:

[0142] The system uses the `score_genes` function of `scanpy` to calculate the functional module score corresponding to each cell. Subsequently, the system sequentially performs ScalePert-Cell virtual cell perturbation calculation, Wasserstein-1 distribution reversal score calculation, ScalePert-Organoid intercellular communication network propagation calculation, and candidate drug tissue-level ranking calculation.

[0143] Appendix Figure 1 The overall architecture of the virtual disturbance analysis system of this invention is shown in the appendix. Figure 1 As can be seen, this system comprises two levels: a virtual cell perturbation module and a virtual organoid perturbation module. The virtual cell perturbation module is used to generate virtual cell states after perturbation of the target gene at single-cell resolution, while the virtual organoid perturbation module is used to propagate the direct perturbation score at the cell type level along the intercellular communication network, thereby forming a tissue-level inverse score.

[0144] This embodiment illustrates that the present invention is not a standalone gene expression scoring method, but rather a continuous computational system encompassing single-cell transcriptome data input, cell state space construction, virtual gene perturbation, module distribution reversal evaluation, intercellular communication propagation, and drug ranking output. This system does not rely on scATAC-seq data or Perturb-seq supervised training data; instead, it directly predicts the target gene perturbation effect based on observational single-cell / mononuclear RNA sequencing data. Therefore, the present invention is applicable to data scenarios such as clinical biopsy samples where real gene perturbation experiments are difficult to conduct.

[0145] Example 2: Single-gene virtual perturbation analysis based on PCA spatially low-expression K-nearest neighbors This embodiment illustrates the calculation process of the single-gene virtual perturbation displacement vector in this invention, and can be referred to in conjunction with the appendix. Figure 2 Appendix Figure 3 and attached Figure 4 Please provide an explanation.

[0146] In this embodiment, the system performs single-gene virtual exhaustion analysis for each of the nine hub genes. The nine hub genes are TLR2, ITGAM, CCL5, HIF1A, VCAM1, JUN, AIM2, IFI16, and PYCARD.

[0147] For the target gene to be perturbed The system extracts the expression vector of a cell from the standardized expression data in all cells.

[0148] The system is the first in the PCA cell state space. A K-nearest neighbor set is constructed from each cell, denoted as _____. .

[0149] In this embodiment, the nearest neighbor number is used for perturbation direction calculation. Nearest neighbor search uses cosine distance.

[0150] After obtaining the K-nearest neighbor set, the system determines the target gene. The expression value is used to filter the low-expression nearest neighbors in the K-nearest neighbor set: ; In other words, The cells in the passage simultaneously satisfy two conditions: first, they belong to the first... The K-nearest neighbor set of each cell in the PCA space; secondly, the target gene. The expression value is lower than the first Each cell.

[0151] when When the set is not empty, the system uses the average direction of the cells in the set relative to the current cell as the local movement direction after the target gene expression decreases, and calculates the target gene. For the first Virtual perturbation displacement vector generated by each cell: ; Based on the perturbation displacement vector, the system calculates the target gene. After the disturbance The virtual position of each cell in PCA space: ; when When the value is empty, the system marks the cell as an invalid perturbation cell under the target gene and sets its perturbation displacement vector accordingly. Zero vectors are used as placeholders. Invalid perturbation cells do not participate in the construction of the pre-perturbation module score distribution and post-perturbation module score distribution corresponding to the target gene. If the number of valid perturbation cells in a certain cell type is lower than the preset minimum number, the system marks the reverse score of the target gene in that cell type as a missing value, or fills it with zeros according to preset rules.

[0152] Virtual location after perturbation Then, the system estimates the functional module score corresponding to the perturbed position using K-nearest neighbor interpolation. In this embodiment, the perturbed module score estimation adopts... The K-nearest neighbor interpolation is performed, using Euclidean distance as the distance metric. Specifically, the system searches for neighbors in the original PCA space. The 10 nearest primitive cells were identified, and functional module scores were estimated based on these neighboring primitive cells. The corresponding module score after the perturbation.

[0153] For any functional module, the system constructs a pre-perturbation module score distribution from the module scores of the effective cells before the perturbation, and constructs a post-perturbation module score distribution from the module scores corresponding to the virtual cell states after the perturbation. The system uses Wasserstein-1 distance to calculate the non-negative displacement between the pre-perturbation and post-perturbation module score distributions. Subsequently, based on the direction of change of the post-perturbation module score distribution relative to the pre-perturbation module score distribution, the system assigns a directional sign to this non-negative displacement, forming a signed inverse score.

[0154] In this embodiment, the reduction of module score is taken as the preset treatment direction, and the treatment direction score is recorded as a negative value; when it is used for candidate target or candidate drug ranking in the future, the system converts the treatment direction score into positive therapeutic reversal intensity.

[0155] To evaluate the statistical support of the reversal score, the system performed 1000 permutations on the correspondence between module scores and cells to construct a random null distribution, and compared the true reversal score with the random null distribution. Simultaneously, the system calculated the Cohen's d effect size; for multiple statistical test results generated from multi-gene, multi-module, and multi-cell types, the Benjamini-Hochberg method was used for FDR correction.

[0156] Appendix Figure 4 The results of single-gene virtual perturbation scoring for nine hub genes in different cell types are presented. According to the calculations in this embodiment, IFI16 depletion produces the strongest global reversal effect in almost all cell types.

[0157] To verify the nearest neighbor number The system's impact on single-gene perturbation results Under different conditions, the virtual perturbation analysis was repeated, and different calculations were performed. Spearman rank correlation coefficients between target gene sequences at different values. (Appendix) Figure 2 Display the number of nearest neighbors The sensitivity analysis results showed that all pairwise Spearman correlation coefficients were greater than 0.9, indicating that changes in the number of nearest neighbors did not significantly alter the target gene sequence.

[0158] To verify the perturbation intensity scaling factor The impact on the results, the system in Under different conditions, the virtual perturbation analysis was repeated, and different calculations were performed. Spearman rank correlation coefficients between target gene sequences at different values. (Appendix) Figure 3 Shows the disturbance intensity scaling factor The sensitivity analysis results showed that all pairwise Spearman correlation coefficients were greater than 0.9, indicating that the single-gene virtual perturbation ordination results of this invention have strong stability with the selection of perturbation intensity scaling factors.

[0159] Example 3: Combined perturbation analysis based on the average superposition of multiple gene perturbations and displacements This embodiment illustrates the calculation process of multi-gene combination perturbation in this invention, and can be combined with the appendix. Figure 5 Please provide an explanation.

[0160] After completing the single-gene virtual perturbation analysis, the system further analyzes combined perturbation objects consisting of two or more target genes. Let the target gene combination corresponding to the combined perturbation object be: ; in, This indicates the number of component genes in the combination, and .

[0161] For the For each cell, the system has calculated the target gene combination according to Example 2. The virtual perturbation shift vectors of each component gene in the data are as follows: ; The system reuses the single-gene virtual perturbation displacement vectors corresponding to each component gene in the same PCA cell state space, and calculates the target gene combination using an equal-weighted averaging method. For the first The combined perturbation displacement vector generated by each cell: ; When a certain component gene In the When there are no effectively low-expression neighbors in a cell, the system will virtually perturb the shift vector of the single gene corresponding to the component gene. This is denoted as the zero vector, and it is included in the calculation of the combined disturbance displacement vector.

[0162] After obtaining the combined disturbance displacement vector Then, the system calculates the target gene combination. After the disturbance The virtual position of each cell in PCA space: ; The virtual position after obtaining the combined perturbation Then, the system estimates the perturbed functional module scores in the same way as the single-gene perturbation.

[0163] Specifically, the system searches for and in the original PCA space The nearest preset number of primitive cells, and the functional module score estimate based on the primitive cells pre-calculated. The corresponding module score after the perturbation.

[0164] For any functional module, the system constructs a pre-perturbation module score distribution from the module scores before the perturbation and a post-perturbation module score distribution from the module scores after the perturbation. The system then calculates the non-negative distributed displacement between the two using the Wasserstein-1 distance. Subsequently, based on the direction of change of the post-perturbation module score distribution relative to the pre-perturbation module score distribution, the system assigns a directional sign to this non-negative distributed displacement, thus obtaining the signed combined reversal score of the target gene combination on the functional module.

[0165] In this embodiment, the reduction of module score is taken as the preset treatment direction, and the treatment direction score is recorded as a negative value; when it is used for subsequent candidate target combination ranking or candidate drug ranking, the system converts the treatment direction score into positive therapeutic reversal intensity.

[0166] When calculating the combined disturbance reversal score, the system is based on the combined disturbance displacement vector. The system selects cells with effective perturbation combinations based on their amplitude. Cells whose perturbation displacement amplitude does not meet the preset conditions are not included in the construction of the module score distribution before and after perturbation. If the number of effective perturbation cells in a certain cell type is lower than the preset minimum number, the system marks the perturbation reversal score of the target gene combination in that cell type as a missing value, or fills it with zero according to preset rules.

[0167] Appendix Figure 5 The results of reversal scoring for bigenic and trigenic combination perturbations in different cell types are presented. (According to the appendix...) Figure 5 The results show that the IFI16+AIM2 combination achieved a reversal score exceeding the sum of individual effects in Myofib cell types, demonstrating a significant synergistic reversal effect. This indicates that the present invention can evaluate not only the perturbation effect following the depletion of a single target gene, but also the combined effect of simultaneous perturbation of multiple target genes.

[0168] Example 4: Tissue-level perturbation propagation based on ligand-receptor communication networks and thermonuclear diffusion This embodiment illustrates the implementation process of the virtual organoid perturbation module of the present invention, and can be referred to in conjunction with the appendix. Figure 6 Appendix Figure 7 Appendix Figure 8 and attached Figure 9 Please provide an explanation.

[0169] In this embodiment, the system constructs an intercellular communication network based on the 11 kidney cell types obtained in Embodiment 1. The set of ligand-receptor relationships used to construct the intercellular communication network includes the following 18 pairs of ligand-receptor relationships:

[0170] The system groups the expression matrix by cell type based on cell type annotation information. For any sent cell type... and receiving cell types The system calculates the sending cell type for each ligand-receptor pair. Mean normalized expression values ​​of ligand genes and percentage of positive cells, as well as recipient cell types. The average normalized expression value of the receptor gene and the proportion of positive cells.

[0171] The system multiplies the four expression parameters mentioned above to obtain the communication score of the ligand-receptor pair in the direction from the sending cell type to the receiving cell type. The communication scores of all ligand-receptor pairs are then summed or weighted to obtain the sending cell type. To receive cell type The communication strength.

[0172] The system further screened statistically significant communication relationships through 200 permutation tests and retained those that met the criteria. The communication relationships between cell types are identified. These filtered communication relationships are then used to construct a directed adjacency matrix. ; In one specific implementation, the system uses a communication adjacency matrix. Normalization is performed. To ensure that subsequent thermonuclear diffusion has a clear propagation meaning in the cell-type network, the system will use a directed communication adjacency matrix. Converted into a nonnegative propagation matrix for thermonuclear diffusion Propagation matrix It can be by It can be obtained through normalization, or it can be obtained from... It is obtained by symmetrizing its transpose matrix.

[0173] Preferably, the system obtains the undirected weighted propagation matrix in the following manner: ; The resulting undirected weighted propagation matrix It possesses numerical symmetry. The system directly uses the undirected weighted propagation matrix. Construct the degree matrix Degree matrix Let be a diagonal matrix, and its first... The diagonal elements are the propagation matrix. The Middle The sum of connection weights corresponding to each cell type, i.e.: ; Based on the propagation matrix The sum-degree matrix DDD, the system constructs the graph Laplacian matrix: ; Thulaplatz matrix Used to represent the topology of cell-type communication networks and as the fundamental matrix for thermonuclear diffusion propagation.

[0174] The system denotes the cell type level direct perturbation scoring matrix as follows: In this embodiment, thermonuclear diffusion parameters The default value is 0.3. The system calculates the thermonuclear diffusion matrix: ; Subsequently, the system obtains the propagated perturbation score matrix using the following formula: ; Thermonuclear diffusion propagates perturbation scores along the cell type dimension, thus resulting in a post-propagation matrix. The cells are still grouped by cell type and the target perturbation object by column.

[0175] It should be noted that, It is not the final tissue-level reversal score, but rather a cell type-level score matrix propagated through intercellular communication networks. The system further reads... The scores of the same target perturbation object after propagation on various cell types are converted into positive reversal intensity for those that conform to the treatment direction. These scores are then aggregated in a way that is averaged, summed, or weighted averaged based on cell type proportions to form a tissue-level reversal score.

[0176] Appendix Figure 7 The communication networks among 11 cell types after significance screening are shown. (See attached image.) Figure 8 The tissue-level results of single-gene perturbations after thermonuclear diffusion are presented. Based on the calculations in this embodiment, the tissue-level single-gene perturbation ranking is as follows: ; Appendix Figure 9 This study presents tissue-level results following the propagation of combined gene perturbations via thermonuclear diffusion.

[0177] Parameter sensitivity analysis: To verify the thermonuclear diffusion parameters The system repeatedly performs thermonuclear diffusion propagation and tissue level score calculation under multiple preset thermonuclear diffusion parameters to influence the organization level ranking.

[0178] In this embodiment, the system respectively: ; Under different conditions, repeat the organization-level propagation analysis and compare different Spearman rank correlation coefficient between the ordering of target perturbation objects under the given value.

[0179] Appendix Figure 6 Thermonuclear diffusion parameters were displayed. The results of the sensitivity analysis showed that different The tissue-level ranking is highly consistent, and all pairwise Spearman correlation coefficients are greater than 0.9, indicating that the tissue-level propagation results of this invention have strong stability for thermonuclear diffusion parameters.

[0180] Example 5: Virtual screening of candidate drugs based on tissue-level reversal scores This embodiment illustrates the application of the present invention in virtual screening of candidate drugs, and can be combined with the appendix. Figure 10 Please provide an explanation.

[0181] In this embodiment, the system receives information on 25 candidate drugs and their known targets. The candidate drug target information includes the candidate drug name and one or more target genes corresponding to that candidate drug. The system matches the target genes of the candidate drugs with the target gene set and determines the target perturbation object corresponding to the candidate drug based on the matching results.

[0182] When a candidate drug has a matching target gene among its target genes, the system treats the candidate drug as a single-gene perturbation object and calls the single-gene virtual perturbation shift, cell type level direct perturbation score, and tissue level reversal score corresponding to the target gene.

[0183] When there are two or more matching target genes in the target gene of a candidate drug, the system treats the candidate drug as a combined perturbation object and calculates the perturbation displacement corresponding to the candidate drug according to the combined perturbation method.

[0184] When the target genes of a candidate drug fail to match any genes in the target gene set or expression matrix, the system marks the candidate drug as an unscored candidate drug and excludes it from the candidate drug ranking, or marks the candidate drug as a missing value for tissue-level reverse scoring.

[0185] In this invention, drug target mapping does not directly predict the binding relationship between drugs and target molecules. Instead, it uses known drug target information to associate candidate drugs with the aforementioned virtual perturbation process. Therefore, drug evaluation can reuse cell state changes, functional module reversal scores, and tissue-level propagation results obtained from single-gene perturbations and combined perturbations.

[0186] Let the first The set of matchable target genes corresponding to each candidate drug is: ; when When the candidate drug contains two or more target genes, the system calculates the perturbation shift corresponding to the candidate drug according to the combined perturbation method. For the first... The cell, the first The perturbation displacement vectors corresponding to each candidate drug are: ; when A specific target gene In the When there are no effectively low-expression neighbors in a single cell, the system will perturb the single-gene shift vector corresponding to the target gene. The zero vector is used in the above calculations, or it is removed from the effective combination components according to the preset rules.

[0187] Preferably, the system uses a zero-vector placeholder method.

[0188] After the candidate drug was perturbed, the first The virtual positions of each cell in PCA space are: ; when When the candidate drug contains only one matching target gene, the system directly calls the single-gene perturbation shift and reversal score results corresponding to that target gene. In this case, the virtual perturbation corresponding to the candidate drug is equivalent to the virtual perturbation of that single target gene.

[0189] when When it is an empty set, it represents the first set. If a candidate drug does not have a matching target gene that can enter the current virtual perturbation process, the system will mark the candidate drug as an unscored candidate drug and exclude it from the candidate drug ranking, or mark its cell type level direct perturbation score and tissue level inverse score as missing values.

[0190] Virtual cell locations after obtaining candidate drug perturbation Subsequently, the system obtains the perturbation module score according to the aforementioned functional module score estimation method. Then, based on cell type annotation, the system divides the virtual cell state of the candidate drug after perturbation according to cell type, constructs the functional module score distribution before and after perturbation within each cell type, and calculates the direct perturbation score of the candidate drug at each cell type level based on Wasserstein-1 distance.

[0191] For each candidate drug, the system inputs its cell type-level direct perturbation score into the intercellular communication network of Example 4. Specifically, the system uses the cell type-level direct perturbation score corresponding to the candidate drug as a direct perturbation score matrix. The column corresponding to the target perturbation object is then processed through the thermal core diffusion matrix. To spread: ; System read The post-spreading score of the corresponding candidate drug is used to convert the score that matches the treatment direction into the positive therapeutic reversal intensity, and then aggregated according to the average, summation or weighted average based on cell type ratio to form the tissue-level reversal score of the candidate drug.

[0192] After calculating the tissue-level reversal scores for all candidate drugs, the system ranks the candidates from highest to lowest based on the intensity of their positive therapeutic reversal. The ranking results are used to determine the priority of candidate drugs within the tissue context corresponding to the current single-cell transcriptome data. When outputting the ranking results, the system retains the drug name, corresponding target gene, cell type-level score, and tissue-level score, allowing users to simultaneously view the overall drug ranking and cell type origin.

[0193] In this embodiment, the candidate drug target mapping includes 25 candidate drugs and their corresponding targets. The system maps each candidate drug to one or more hub gene perturbation objects and calculates its cell type-level direct perturbation score and tissue-level reversal score, respectively. (Appendix) Figure 10 The results of tissue-level reversal scoring and ranking of 25 candidate drugs after thermonuclear diffusion are presented. According to the calculation results of this embodiment, VX-765 ranks first among the 25 candidate drugs.

[0194] From the appendix Figure 10 As can be seen, VX-765 ranked first among 25 candidate drugs, indicating that perturbation of its corresponding target can produce a strong tissue-level therapeutic reversal effect in the context of cell communication networks. This result supports the practical application value of this invention for virtual screening and prioritization of candidate drugs.

[0195] It should be noted that some English labels in the attached figures are commonly used technical terms, gene names, drug names, cell type abbreviations, or algorithm output fields. For ease of understanding, their meanings are explained as follows in the manual: Cell Type or cell_type indicates cell type; target indicates the target perturbation object; Hub Gene indicates the candidate hub gene; Wasserstein Reversal indicates the reversal score calculated based on Wasserstein-1 distance; TissueReversal Score indicates the tissue-level reversal score obtained by propagation and aggregation through intercellular communication networks; Reversalby Cell Type indicates the reversal score displayed by cell type; Permutation-filtered indicates the results after permutation test screening; Spearman ρ indicates the Spearman rank correlation coefficient; values>0.9 indicate robustness, meaning that a correlation coefficient greater than 0.9 indicates that the sorting results are robust.

[0196] The following is a table of English and Chinese labels for the attached diagram:

[0197] The following is an appendix Figure 4 - Appendix Figure 10 Cell type reference table used for the horizontal or vertical axis:

Claims

1. A target virtual perturbation analysis system based on single-cell transcriptome data, characterized in that, include: The data input module is used to obtain the single-cell transcriptome expression matrix and the target gene set; The space construction module is used to construct a cell state space based on the single-cell transcriptome expression matrix; The virtual cell perturbation layer is used to perform virtual perturbation analysis on the target gene in the cell state space. It calculates the virtual perturbation displacement vector based on the neighborhood state of each cell in the cell state space to obtain the perturbated virtual cell state. The perturbation evaluation layer is used to calculate the perturbation evaluation score of the target gene based on the changes in the virtual cell state compared to before the perturbation. The result output module is used to output the evaluation result for the target gene based on the perturbation evaluation score.

2. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 1, characterized in that, The space construction module is specifically used for constructing the cell state space as follows: Cells with the number of detected genes within a preset range and the proportion of mitochondrial-derived reads below a preset threshold are retained, while genes that are underexpressed in the cell population are removed. The total expression level of each cell in the quality-controlled single-cell transcriptome expression matrix was normalized and logarithmically transformed, and a predetermined number of hypervariable genes were screened from the standardized expression matrix. The expression values ​​of the hypervariable genes are scaled, centered, and variance standardized, and principal component analysis is performed to reduce dimensionality, extracting the first few principal components as the cell state space.

3. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 1, characterized in that, The virtual cell perturbation layer is specifically used to calculate the virtual perturbation displacement vector: The expression vector of the target gene in all cells is extracted from the standardized expression data obtained by standardizing the single-cell transcriptome expression matrix. In the cell state space, a K nearest neighbor set with a preset number of nearest neighbors is constructed for any current cell, and nearest neighbor cells with target gene expression values ​​lower than those of the current cell are selected from the K nearest neighbor set to form a low-expression nearest neighbor set; The coordinates of the current cell are subtracted from the coordinates of each neighboring cell in the low-expression nearest neighbor set in the cell state space to obtain the direction vector. The direction vector is then averaged with equal weights and multiplied by the perturbation intensity scaling factor to calculate the virtual perturbation displacement vector generated by the target gene on the current cell. The original coordinates of the current cell in the cell state space are added to the virtual perturbation displacement vector to obtain the position after virtual perturbation, thus forming the virtual cell state.

4. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 3, characterized in that, The virtual cell perturbation layer is also used to handle invalid cells at the boundary: When the set of low-expression nearest neighbors corresponding to the current cell is empty, the current cell is marked as an invalid perturbation cell under the target gene, and the virtual perturbation displacement vector corresponding to the current cell is occupied by a zero vector; The invalid perturbation cells are set not to participate in the construction of the pre-perturbation module score distribution and post-perturbation module score distribution required for subsequent perturbation evaluation; If the number of effective perturbed cells in a certain cell type is lower than the preset minimum number, the target gene is marked as a missing value or filled with zero according to the preset rules as the reversal score of the perturbed evaluation score in that cell type.

5. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 3, characterized in that, The target gene set contains two or more target genes as a combination of target genes. The virtual cell perturbation layer, when calculating the combined perturbation, is specifically used for: The virtual perturbation displacement vectors corresponding to each target gene are reused in the same cell state space; For the current cell, the average value of the virtual perturbation displacement vectors corresponding to each component gene in the target gene combination is calculated to obtain the combined perturbation displacement vector; If a component gene has no effectively low-expression neighbors in the current cell, the virtual perturbation displacement vector corresponding to the component gene is recorded as the zero vector and used in the combination calculation. The original coordinates of the current cell in the cell state space are added to the combined perturbation displacement vector to obtain the virtual position after combined perturbation.

6. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 1, characterized in that, The data input module is also used to acquire cell type annotation information and functional module gene sets; The perturbation evaluation layer is specifically a virtual organoid perturbation layer, which is used to divide the virtual cell state into communities according to cell type based on the cell type annotation information, and to construct a pre-perturbation functional module score distribution composed of cells that meet the preset effective perturbation conditions, and a post-perturbation functional module score distribution composed of the corresponding virtual cell state for each cell type. The Wasserstein-1 distance is used to calculate the distribution displacement between the functional module score distributions before and after the perturbation. Based on whether the change direction of the distribution after the perturbation relative to the distribution before the perturbation conforms to the preset treatment direction, the distribution displacement is assigned a treatment direction symbol or a non-treatment direction symbol to generate a signed inverse score, which serves as a direct perturbation score at the cell type level.

7. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 6, characterized in that, The functional module scores of the original cells involved in the pre-perturbation functional module score distribution and the post-perturbation functional module score distribution are calculated in advance by the spatial construction module or the data input module based on the functional module gene set. The algorithm used includes at least one of GSVA, AUCell algorithm, UCell algorithm or the AddModuleScore method of the Seurat platform.

8. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 6, characterized in that, The data input module is also used to obtain a set of ligand-receptor relationships; The virtual organoid perturbation layer is specifically used in constructing intercellular communication networks for: The average normalized expression values ​​of ligand genes and receptor genes and the proportion of positive cells are extracted according to cell type. The communication score of each ligand-receptor pair is calculated by multiplying the average expression value of ligands and the proportion of positive ligands in the sending cell type with the average expression value of receptors and the proportion of positive receptors in the receiving cell type. The directed communication adjacency matrix is ​​obtained by summing or weighted summing the communication scores of all ligand-receptor pairs. The directed communication adjacency matrix and its transpose matrix are symmetricized to construct an undirected weighted propagation matrix. The sum of connection weights for each cell type is calculated based on the undirected weighted propagation matrix to construct a degree matrix. The degree matrix is ​​a diagonal matrix. The graph Laplacian matrix is ​​constructed by subtracting the undirected weighted propagation matrix from the degree matrix.

9. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 8, characterized in that, The virtual organoid perturbation layer, when propagating the direct perturbation score, is specifically used for: The signed inverse scores of multiple functional modules under each cell type are combined into a comprehensive direct perturbation score by means of equal weighted average, preset weighted average, average of scores of only the treatment direction module, or weighted summation according to disease relevance weight, so as to construct a cell type level direct perturbation score matrix. A thermonuclear diffusion matrix is ​​constructed using the aforementioned Graph Laplacian matrix and preset thermonuclear diffusion parameters. The propagated perturbation score matrix is ​​obtained by multiplying the thermonuclear diffusion matrix with the cell type level direct perturbation score matrix. The scores of the target perturbation object in the propagated perturbation score matrix are read, the scores that conform to the treatment direction are converted into positive reversal intensity, and the tissue-level reversal scores are aggregated in the form of averaging, summation or weighted averaging based on cell type proportion.

10. The target virtual perturbation analysis system based on single-cell transcriptome data according to claim 9, characterized in that, The data input module is also used to acquire candidate drug target information; When the result output module sorts candidate drugs based on the candidate drug target information, it is specifically used for: Match the target genes of the candidate drugs with the set of target genes; When there is a matching target gene among the target genes, the candidate drug is used as a single gene perturbation object and its corresponding tissue-level reversal score is invoked. When there are two or more matching target genes in the target gene, the candidate drug is used as a combined perturbation object and the corresponding tissue-level reversal score is invoked. When none of the target genes match the genes in the target gene set or expression matrix, the candidate drug is marked as an unscored candidate drug and excluded from the sorting, or the tissue-level reverse score of the candidate drug is marked as a missing value. The tissue-level reversal scores of each candidate drug that conform to the treatment direction are converted into positive therapeutic reversal intensity, and the candidate drugs are sorted from high to low according to the positive therapeutic reversal intensity as the candidate drug ranking.

Images