A Method and System for Drug Response Cell Population Ranking Based on Single-Cell Transcriptome Data

Abstract

This invention provides a method and system for ranking drug response cell populations based on single-cell transcriptome data, relating to the field of drug response cell population ranking. The method includes: constructing a target gene regulatory network (GRN) for each cell population based on scRNA-seq data under disease states of known cell population types; virtually knocking out target points in the GRN based on drug target information to construct a target gene regulatory network (tpGRN) for each cell population after virtual target knockout; training a low-dimensional representation of each network node in the GRN and the tpGRN using manifold alignment; calculating the distance of each network node in the GRN and the tpGRN using Euclidean distance; comprehensively considering the changing trends of drug targets, 2-hop nodes, and edges of 2-hop nodes in the GRN and tpGRN to score the drug response of different cell populations, generating drug perturbation scores to rank the drug responses of cell populations and determine the ranking results. This invention can improve the accuracy of inferring drug response cell populations.

Description

Technical Field

[0001] This invention relates to the field of drug response cell population sequencing, and in particular to a method and system for drug response cell population sequencing based on single-cell transcriptome data. Background Technology

[0002] Single-cell RNA-seq (scRNA-seq) technology enables the determination of all gene expression in each single cell and has been widely used to study cellular heterogeneity. Given the high heterogeneity of cellular responses to drugs, identifying the most drug-sensitive cell populations using scRNA-seq data is crucial for understanding drug mechanisms of action and guiding drug development, and has attracted considerable attention. Currently, the most common approach is to score and rank drug-responsive cell populations based on the number of differentially expressed genes (DEGs) between paired cell populations under different conditions. Recently, some researchers have also proposed a method called Augur for perturbation scoring and ranking of cell populations. This method utilizes the separability of the same cell population under different conditions in a high-dimensional space to prioritize cell types most sensitive to biological perturbations. However, drugs typically rely on binding to target proteins to exert their therapeutic effects. Clearly, these methods do not consider prior knowledge of drug targets, leading to inaccurate inferences about drug-responsive cell populations. Furthermore, existing strategies and methods require scRNA-seq data from two different states, severely limiting their application in scRNA-seq data with only disease states available. To date, there is no method for predicting drug-responsive cell populations using only scRNA-seq data from disease states. Summary of the Invention

[0003] The purpose of this invention is to provide a method for ranking drug response cell populations based on single-cell transcriptome data, in order to solve the problem of low accuracy in inferring drug response cell populations.

[0004] To achieve the above objectives, the present invention provides the following solution:

[0005] A method for ranking drug response cell populations based on single-cell transcriptome data, comprising:

[0006] Obtain scRNA-seq data of disease states and drug target information for known cell population types;

[0007] Based on the scRNA-seq data, a target gene regulatory network (GRN) for each cell population was constructed using gene co-expression relationships.

[0008] Based on the drug target information, the target in the GRN is virtually knocked out to construct the target gene regulatory network tpGRN for each cell population after virtual knockout of the target.

[0009] The low-dimensional representation of each network node in the GRN and the tpGRN is trained by manifold alignment; the network node is a gene node.

[0010] Based on the low-dimensional representation, the distance of each network node in the GRN and the tpGRN is calculated using Euclidean distance;

[0011] Based on the distance, the drug response of different cell populations is scored by comprehensively considering the changing trends of the drug target, 2-hop node and edge of the GRN and the tpGRN, and a drug perturbation score is generated.

[0012] The cell populations are ranked according to the drug perturbation score to determine the ranking of the drug-response cell populations; the ranking of the drug-response cell populations is used to characterize the degree of response of the cell populations to drug perturbation.

[0013] Optionally, based on the scRNA-seq data, a target gene regulatory network (GRN) for each cell population is constructed using gene co-expression relationships, specifically including:

[0014] Construct a scRNA-seq expression matrix based on the scRNA-seq data;

[0015] δc cells are randomly selected from each cell type in the scRNA-seq expression matrix, and the gene expression matrix corresponding to the selected δc cells is extracted, and this process is repeated S times; where δ is the selection ratio and c is the number of cells.

[0016] Feature selection is performed on the genes in the gene expression matrix to construct a gene set of GRNs; the gene set includes the top 2000 highly variable genes in the scRNA-seq data, genes corresponding to transcription factors in the AnimalTFDB database, and drug target genes in the DGIdb database;

[0017] Principal component analysis was used to determine the low-dimensional representation of each gene in the gene set and to estimate the regression coefficients between genes as the weights of the edges in the GRN.

[0018] Define any gene expression matrix, and use principal component analysis to regress the target gene on multiple potential covariates to construct the relationship between the target gene and the regulatory gene;

[0019] A gene adjacency matrix is ​​constructed based on the relationship between the target gene and the regulatory gene and the weights thereon.

[0020] Tensor component analysis was used to combine the adjacency matrices of multiple genes into a third-order tensor χ of S×T×R; where S is the number of GRNs, T is the number of target genes, and R is the number of regulatory genes.

[0021] Extract the main feature patterns shared by multiple GRNs, decompose and approximate the third-order tensor χ, and generate a new third-order tensor.

[0022] Based on the new third-order tensor, each weight is normalized by dividing it by the largest absolute value among all weight absolute values ​​to construct the final adjacency matrix for each cell population; the final adjacency matrix is ​​the target gene regulatory network GRN.

[0023] Optionally, the relationship between the target gene and the regulatory gene is as follows:

[0024] Y = X′β′ + ε;

[0025] Where Y is the expression vector of the target gene; X′ is the expression vector of the regulatory genes corresponding to the top-ranking principal components; β′ is the regression coefficient of the latent covariates after transformation; and ε is the random error.

[0026] Optionally, extract the main feature patterns shared by multiple GRNs, decompose and approximate the third-order tensor χ, and generate a new third-order tensor, specifically including:

[0027] Extract the main feature patterns shared by multiple GRNs, and use the formula The third-order tensor χ is decomposed and approximated; where s r Let r be the vector of factors r in the number dimension of GRN; For outer product operation; t r Let r be the vector of factor r in the target gene dimension; r The vector that regulates factor r in the gene dimension;

[0028] The first five factors are selected to characterize the key components of the third-order tensor χ, and a new third-order tensor is generated by averaging the accumulated components.

[0029] Optionally, based on the drug target information, the target in the GRN is virtually knocked out to construct the target gene regulatory network tpGRN for each cell population after virtual target knockout, specifically including:

[0030] Extract the row containing the drug target gene from the final adjacency matrix and set the corresponding value of the row containing the drug target gene to 0 to virtually construct the drug-perturbed adjacency matrix; the drug-perturbed adjacency matrix is ​​the tpGRN.

[0031] Optionally, a low-dimensional representation of each network node in the GRN is trained in the GRN and the tpGRN through manifold alignment, specifically including:

[0032] The final adjacency matrix and the drug-perturbed adjacency matrix are combined into a joint matrix;

[0033] Using manifold alignment, the joint matrix is ​​projected onto a shared low-dimensional latent space to determine the low-dimensional representation of each gene node in the GRN and the tpGRN, and the Euclidean distance between each matched gene node is interpreted as the perturbation effect of the drug on the target cell type.

[0034] Optionally, the joint matrix is:

[0035]

[0036] Where W is the joint matrix; A′ is the final adjacency matrix; A′ PB λ is the adjacency matrix after drug perturbation; λ is the adjustment parameter; I is the identity matrix describing the correspondence between gene nodes in the GRN and tpGRN.

[0037] Optionally, the drug perturbation score is:

[0038]

[0039] Among them, D Tar W represents the distance between the target gene and the GRN and the tpGRN. out The weight of the edge leaving the gene node in the GRN; Deg Tar The degree of the target gene in the GRN; N is the total number of genes regulated by the target gene; n is the gene number; D n W represents the distance between the downstream gene and the GRN and the tpGRN. in The weight of the edge entering the gene node in the GRN; Deg n The degree of the downstream gene in the GRN.

[0040] A drug response cell population ranking system based on single-cell transcriptome data, comprising:

[0041] The data and information acquisition module is used to acquire scRNA-seq data and drug target information in disease states with known cell population types.

[0042] The GRN construction module is used to construct the target gene regulatory network (GRN) for each cell population based on the scRNA-seq data and the co-expression relationship of genes.

[0043] The tpGRN construction module is used to virtually knock out the target in the GRN according to the drug target information, and construct the target gene regulatory network tpGRN after virtual knockout of the target for each cell population.

[0044] A manifold alignment module is used to train a low-dimensional representation of each network node in the GRN in the GRN and the tpGRN through manifold alignment; the network node is a gene node.

[0045] The Euclidean distance calculation module is used to calculate the distance of each network node in the GRN and the tpGRN based on the low-dimensional representation using Euclidean distance.

[0046] The drug perturbation score generation module is used to score the drug response of different cell populations based on the distance, taking into account the changing trends of the drug target, 2-hop node and the edge of the 2-hop node in the GRN and the tpGRN, and to generate a drug perturbation score.

[0047] The sorting module is used to sort the drug response of the cell population according to the drug perturbation score and determine the drug response cell population sorting result; the drug response cell population sorting result is used to characterize the degree of response of the cell population to drug perturbation.

[0048] According to specific embodiments provided by the present invention, the following technical effects are disclosed: The present invention provides a method and system for ranking drug response cell populations based on single-cell transcriptome data. It constructs a target gene regulatory network (GRN) for each cell population using only scRNA-seq data under disease conditions, and virtually knocks out target sites in the GRN using drug target information, constructing a target gene regulatory network (tpGRN) for each cell population after virtual target site knockout. Based on the GRN and tpGRN, the drug response of different cell populations is scored, and finally, the drug response of cell populations is ranked according to the drug perturbation score to determine the ranking result of drug response cell populations. The present invention improves the accuracy of inferring drug response cell populations based on prior knowledge of drug targets, and can obtain the drug response of cell populations using only scRNA-seq data under disease conditions, expanding the application scope of scRNA-seq data only under disease conditions. Attached Figure Description

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

[0050] Figure 1 The flowchart of the drug response cell population ranking method based on single-cell transcriptome data provided by the present invention is shown below.

[0051] Figure 2 This is a diagram of the scRank working architecture provided by the present invention; wherein, Figure 2 (a) in the diagram is a schematic of the data input layer; Figure 2 (b) in the diagram is a schematic diagram of the cell type sorting layer; Figure 2 (c) in the diagram is a schematic diagram of the data output layer; Figure 2 (d) in the diagram is a schematic diagram of GRN construction; Figure 2 (e) in the diagram is a schematic diagram of calculating the distance between gene nodes before and after the perturbation; Figure 2 (f) in the diagram is a schematic diagram for calculating the diffusion effect at the target center;

[0052] Figure 3 This is a schematic diagram illustrating the performance test of scRank on simulated data; where, Figure 3 (a) in the diagram is a schematic of SERGIO generating simulation data; Figure 3 (b) in the diagram represents application scenario 1, which is a simulated data diagram showing the orderly decrease in the activity of the co-expression module of the same gene. Figure 3 (c) in the diagram represents application scenario 2, which is a schematic diagram of simulated data showing significant differences in the activity of gene co-expression modules. Figure 3 (d) in the diagram represents application scenario 3, which is a schematic diagram of generating simulated data with cell type-specific gene co-expression module activity.

[0053] Figure 4 This is a schematic diagram illustrating the application of scRank in mouse medulloblastoma data; where, Figure 4 (a) in the figure is a schematic diagram of the prediction results of scRank in mouse medulloblastoma; Figure 4 (b) in the figure is a schematic diagram of cell populations that are sensitive to and resistant to vismodegib, as confirmed by the literature; Figure 4 (c) in the diagram is a schematic diagram of the gene regulatory network of Node_A and Node_B in the disease state; Figure 4 (d) in the figure is a schematic diagram of GRN in two tumor subtypes and their differences before and after drug administration; Figure 4 (e) in the figure is a schematic diagram of the GO enrichment results of gene modules M1 and M2; Figure 4 (f) in the diagram is a schematic diagram of differentially expressed genes in Node_A and Node_B under disease conditions; Figure 4 A schematic diagram of the pathways that are significantly activated in (g)Node_B and Node_A; Figure 4 (h) in the figure represents the prognostic analysis of medulloblastoma patients classified according to the characteristic genes of Node_A and Node_B.

[0054] Figure 5 This is a schematic diagram illustrating the application of scRank in dorsolateral prefrontal cortex data from patients with severe depression; where, Figure 5 (a) shows the scRNA-seq data of patients with severe depression and a schematic diagram of the target SLC6A4 of the therapeutic drug fluoxetine. Figure 5 (b) in the figure is a schematic diagram of the scRank ranking results of all cell populations on the severe depression patient dataset; Figure 5 (c) in the diagram is a schematic diagram of the gene regulatory network and cluster analysis results of Ex_9 and Inhib_5 in the disease state; Figure 5 (d) in the diagram is a schematic diagram of the pathway enrichment results of gene co-expression modules M1, M2, M3, and M4; Figure 5 (e) in the diagram is a schematic diagram of the activity distribution of the M2 gene co-expression module in all cell populations in the MDD; Figure 5 (f) in the figure is a schematic diagram of the differentially expressed genes of Ex_9 in normal and MDD; Figure 5 (g) in the diagram represents a pathway that is significantly enriched in Ex_9; Figure 5 (h) in the figure represents the spatial distribution of different cell populations in the mouse forebrain in MDD; Figure 5 (i) in the figure is a spatial distribution map of different cell populations in the human dorsolateral prefrontal cortex in MDD;

[0055] Figure 6 A comparison graph of scRank with other methods; Figure 6 (a) in the figure is a schematic diagram comparing the performance of scRank, Augur and Wilcox using simulated data; Figure 6 (b) in the figure is a schematic diagram of five sets of scRNA-seq datasets with real drug perturbation; Figure 6 (c) in the figure is a schematic diagram comparing the performance of scRank, Augur, Wilcox, MAST, bimod and DESeq2 using simulated data; Figure 6 (d) in the figure represents the box plot of the prediction results from different methods. Detailed Implementation

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

[0057] The purpose of this invention is to provide a method and system for ranking drug response cell populations based on single-cell transcriptome data, which can improve the accuracy of inferring drug response cell populations.

[0058] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0059] Example 1

[0060] like Figure 1 As shown, this invention provides a method for drug response cell population ranking based on single-cell transcriptome data, including:

[0061] Step 101: Obtain scRNA-seq data of disease status and drug target information for known cell population types.

[0062] Step 102: Based on the scRNA-seq data, construct the target gene regulatory network (GRN) for each cell population using the co-expression relationships of genes.

[0063] In practical applications, step 102 specifically includes: constructing an scRNA-seq expression matrix based on the scRNA-seq data; randomly selecting δc cells from each cell type in the scRNA-seq expression matrix, extracting the gene expression matrix corresponding to the selected δc cells, and repeating this process S times to obtain S GRNs; where δ is the selection ratio and c is the number of cells; performing feature selection on the genes in the gene expression matrix to construct a gene set for the GRNs; the gene set includes the top 2000 highly variable genes in the scRNA-seq data, genes corresponding to transcription factors in the AnimalTFDB database, and drug target genes in the DGIdb database; using principal component analysis to determine the low-dimensional representation of each gene in the gene set and estimating the regression coefficients between genes as the GRNs. The weights of edges in N are defined; for any gene expression matrix, principal component analysis is used to regress the target gene on multiple potential covariates to construct the relationship between the target gene and the regulatory gene; a gene adjacency matrix is ​​constructed based on the relationship between the target gene and the regulatory gene and the weights; multiple gene adjacency matrices are combined into a third-order tensor χ of S×T×R using tensor component analysis; where S is the number of GRNs, T is the number of target genes, and R is the number of regulatory genes; the main feature patterns shared by multiple GRNs are extracted, and the third-order tensor χ is decomposed and approximated to generate a new third-order tensor; based on the new third-order tensor, each weight is divided by the largest absolute value among all weight absolute values ​​for normalization, and the final adjacency matrix for each cell population is constructed; the final adjacency matrix is ​​the target gene regulatory network GRN.

[0064] As an optional embodiment of the present invention, the construction process of the gene regulatory network GRN is as follows:

[0065] For a given cell population i (i∈{1,2,...,C}, the scRNA-seq expression matrix M i [n×c] (n genes and c cells), where C represents the total number of cell types in the dataset. Considering the heterogeneity between cell populations, the drug response cell population ranking method (scRank) based on single-cell transcriptome data provided in this invention first starts from M... i δc cells are randomly selected and randomized t times to simulate a gene expression matrix O = [X1, ..., Xc cells] for δc cells. t [Used to construct cell population-specific GRNs. Where δ represents the selection ratio, δ∈(0,1).]

[0066] To preserve meaningful biological data while reducing computational costs, scRank performed feature selection on the genes in the sequencing data. These featured genes mainly included (1) the top 2000 highly variable genes (HVGs) in the scRNA-seq data, (2) genes corresponding to transcriptional factors (TFs) in the AnimalTFDB database (http: / / bioinfo.life.hust.edu.cn / AnimalTFDB), and (3) target genes of drugs or inhibitors in the DGIdb database (https: / / dgidb.genome.wustl.edu / ). In addition, scRank removed mitochondrial, ribosome, and other redundant genes to preserve intercellular differences and meaningful biological processes. Therefore, the final gene set used to construct GRN was G = [HVG + TF + Tar], with the number represented by p.

[0067] Based on the predefined cells and genes mentioned above, scRank uses principal component analysis (PCA) to obtain a low-dimensional representation of each gene and estimates the regression coefficients between genes as the weights of edges in the GRN. Given a gene expression matrix X[n×p], Y n×1 =(y1, ..., y n ) T X represents the expression vector of the target gene. n×(p-1) = (x1, ... x) n ) T This represents the expression vector of the remaining p-1 regulated genes. Assuming Y depends on X, the relationship between the target gene and other genes can be modeled using the following formula:

[0068] Y = Xβ + ε (1).

[0069] in Let represent the regression coefficient vector in p-1 dimensions in real space, and ε represent the random error. To reduce computational cost and overfitting, scRank uses PCA to regress the target gene on k potential covariates.

[0070] Assume X′ n×k =X n×(p-1) V (p-1)×k X represents the top k principal components selected. n×(p-1) Given the original data, the regression problem can be transformed into:

[0071] Y = X′β′ + ε (2).

[0072] V in the above formula (2) (p-1)×k X is obtained through singular value decomposition and represents a load matrix with p-1 rows and k columns, where each column represents the load of a principal component. This represents the transformed regression coefficients, estimated using ordinary least squares.

[0073] The correlation between the target gene and other genes is calculated using the following formula:

[0074]

[0075] In the above formula (3) and These represent the parameter estimates of the original regression coefficient β and the transformed regression coefficient β′, respectively.

[0076] The first k principal component loads V k As a weight, the regression coefficient β′ of the transformed potential covariate can be transformed into the regression coefficient β of the original variable. Regression through dimensionality reduction can alleviate the multicollinearity problem and increase computational efficiency.

[0077] For each randomly selected feature gene in the iterative regression G, the final result is an adjacency matrix containing t feature gene pairs with their relationships and strengths in cell type i.

[0078] In order to preserve important biological signals in GRN, scRank retains only the top 5% of absolute weight values ​​and removes autocorrelation coefficients to construct the final gene adjacency matrix A.

[0079] To improve the robustness and maintain heterogeneity of GRNs, scRank extracts multiple GRNs from randomly selected data and uses Tensor Component Analysis (TCA) to combine them into a single matrix. Specifically, t GRNs are combined into a third-order tensor χ of size S×T×R, where S, T, and R correspond to the number of GRNs, the number of target genes, and the number of regulatory genes, respectively. To reduce noise in the gene adjacency matrix A, the main shared feature patterns among the multiple GRNs are extracted, and χ is decomposed and approximated using CANDECOMP / PARAFAC, as shown in the following formula:

[0080]

[0081] in, s represents the outer product. r , t r and r r The vector representing factor r contains the loads of each element in each dimension of the tensor. In practical applications, scRank selects the top five factors to represent the key components of χ, obtains a new tensor χ′ by averaging the sums, and normalizes each weight by dividing it by its maximum absolute value to obtain the final GRN for a given cell type i, denoted as A′.

[0082] Step 103: Virtually knock out the target in the GRN according to the drug target information, and construct the target gene regulatory network tpGRN after virtual knockout of the target for each cell population.

[0083] In practical applications, step 103 specifically includes: extracting the row containing the drug target gene from the final adjacency matrix and setting the corresponding value of the row containing the drug target gene to 0, thereby virtually constructing the drug-perturbed adjacency matrix; the drug-perturbed adjacency matrix is ​​the tpGRN.

[0084] Step 104: Train the low-dimensional representation of each network node in the GRN in the GRN and the tpGRN by manifold alignment; the network node is a gene node.

[0085] In practical applications, step 104 specifically includes: combining the final adjacency matrix and the drug-perturbed adjacency matrix into a joint matrix; using manifold alignment, projecting the joint matrix onto a shared low-dimensional latent space, determining the low-dimensional representation of each gene node in the GRN and the tpGRN, and interpreting the Euclidean distance between each matched gene node as the perturbation effect of the drug on the target cell type.

[0086] As an optional embodiment of the present invention, the steps of GRN virtual knockout and manifold alignment are as follows:

[0087] To simulate the effect of drugs or inhibitors on GRNs, scRank extracts the rows containing the drug target genes from the adjacency matrix A′ and sets all values ​​thereto to 0, thereby virtually constructing the drug-perturbed adjacency matrix A′. PB This refers to tpGRN. Manifold alignment is used to jointly project the two high-dimensional adjacency matrices, the original GRN and tpGRN, into a shared low-dimensional latent space, interpreting the distance between each matched gene node as a perturbation effect of the drug on the target cell type.

[0088] Manifold alignment projects two datasets generated by similarity processes into a new low-dimensional space that not only preserves the minimum distance between matching points but also retains their original structure. Therefore, a larger distance between gene nodes indicates a greater change in the topological structure of that node in the two datasets. scRank, through... and The manifold alignment is achieved by combining these into a joint matrix W, as shown in the following formula:

[0089]

[0090] Among them, A′ and A′ PB The variable names in the middle row and column are both gene sets G used to construct the GRN, I represents the identity matrix describing the correspondence between gene nodes in two GRNs, and λ is the adjustment parameter.

[0091] Based on shared low-dimensional projection, scRank uses the Laplacian feature map problem to find d smallest non-zero eigenvalues. Let E 2p×d = (e1, ..., e 2p ) represents the eigenvector, and is expressed by D(p)=||e p -e 2p ||=(d1,...,d p Calculate the Euclidean distance between two gene nodes, where e 2p For the eth 2p There are 10 feature vectors. The Euclidean distance here reflects the impact of target perturbation on gene nodes, simulating the effect of drug action on the gene network.

[0092] Step 105: Based on the low-dimensional representation, calculate the distance of each network node in the GRN and the tpGRN using Euclidean distance.

[0093] Step 106: Based on the distance, the drug response of different cell populations is scored by comprehensively considering the changing trends of the drug target, 2-hop node and edge of the GRN and tpGRN, and a drug perturbation score is generated.

[0094] Step 107: Sort the drug response of the cell population according to the drug perturbation score, and determine the drug response cell population ranking result; the drug response cell population ranking result is used to characterize the degree of response of the cell population to drug perturbation.

[0095] As an optional embodiment of the present invention, the steps of drug-responsive cell population sorting are as follows:

[0096] To more accurately infer the degree of cell population response to drugs, scRank introduces the network diffusion effect. Given a target gene Tar, the influence of drug perturbation propagates from the target gene to downstream genes in the 2-hop of the target gene, and this influence is positively correlated with edge weights and distances between nodes. The specific calculation formula is as follows:

[0097]

[0098] Where D Tar and D n These are obtained through formulas (1)-(5), which represent the distances between the target gene and downstream genes in the GRN and tpGRN before and after virtual knockout, respectively, while Deg represents the degree of the corresponding gene node in the GRN network, and W... in and W out , representing the weights of the edges entering and leaving nodes in the GRN, respectively. Therefore, in the above formulas (1)-(6), the first term represents the normalized drug perturbation effect of the target gene, the second term represents the normalized diffusion effect of the target gene 1-hop nodes, and the third term represents the normalized diffusion effect of the target gene 2-hop nodes. By combining the direct influence of the drug on the gene nodes in step (3) with the diffusion effect of this perturbation in the network, a comprehensive drug perturbation score, or "score," is calculated. Then, this score is used to score and rank the drug response of the cell population.

[0099] The drug response cell population ranking method (scRank) based on single-cell transcriptome data provided in this invention mainly consists of three layers: Figure 2 (a) Data input layer, Figure 2 (b) Cell type sorting layer and Figure 2 The (c) data output layer, in which the cell type sorting layer is further divided into... Figure 2 (d), (e), and (f) in the text.

[0100] The input to scRank consists of scRNA-seq data of disease states for known cell types and drug target information. scRank's cell population ranking algorithm comprises three main parts: First, it constructs a Gene Receiver (GRN) for each cell population based on gene co-expression relationships; second, it virtually knocks out the target in the GRN based on drug target information, constructing a virtual knockout tpGRN for each cell population, and trains the low-dimensional representation of each network node in the GRN before and after virtual knockout using manifold alignment; third, it calculates the distance between each network node before and after virtual knockout using Euclidean distance, comprehensively considering the changes in drug target, 2-hop nodes, and their edges in the GRN before and after virtual knockout to score the drug response of different cell populations. The output of scRank is the rank-sum ranking result, representing the degree of response of the cell population to drug perturbation.

[0101] Example 2

[0102] To implement the method corresponding to Embodiment 1 above, this invention generates scRNA-seq simulated data using a single-cell gene expression simulator (SERGIO) and a predefined GRN module. Based on the scientific assumption that the effect of drug perturbation depends on gene activity, drug-responsive and non-responsive cell populations, i.e., positive and negative samples, are defined. Considering that the perturbation effect of densely connected subnetworks dominates compared to sparsely connected subnetworks, three application scenarios were designed using the simulated data to benchmark scRank, thereby evaluating its predictive performance in identifying drug-responsive cell populations.

[0103] (1) Simulation data construction:

[0104] To construct simulated data consistent with real-world application scenarios, the "SERGIO" Python package and predefined GRN modules were used to automatically generate gene expression profile matrices for specific cell populations. The GRN contains four gene co-expression modules, each containing 25 genes. The co-expression intensity within each module was adjusted by modifying productivity and interaction strength parameters to simulate the expression activity of cell population-specific gene modules. By adjusting the parameters and GRN, three different scRNA-seq simulation datasets and their corresponding three application scenarios were designed to benchmark scRank. Figure 3 As shown, Figure 3 This is a schematic diagram illustrating the performance test of scRank on simulated data. Figure 3(a) in the diagram is a schematic of SERGIO generating simulation data. Figure 3 (b) in the diagram represents application scenario 1, illustrating the simulated data showing the orderly decrease in the activity of co-expression modules of the same gene. Figure 3 (c) in the diagram represents application scenario 2, illustrating simulated data showing significant differences in the activity of gene co-expression modules. Figure 3 (d) in the diagram represents application scenario 3, which is a schematic diagram of generating simulated data with cell type-specific gene co-expression module activity.

[0105] Application Scenario 1 aims to generate simulated data showing an ordered decrease in the activity of co-expression modules of the same gene. The specific steps are as follows: Four datasets with different cell numbers are generated using SERGIO simulation, corresponding to cell populations of 1000, 3000, 5000, and 10000 respectively. For each dataset, four cell populations (A, B, C, D) are simulated, each containing 100 genes, corresponding to four gene modules (M1, M2, M3, M4). These include three co-expression modules with similar activities (M2, M3, and M4) and one co-expression module (M1) whose activity gradually decreases across the four cell populations. Assuming the drug's target gene is located in gene module M1, the actual drug response order across the four cell populations should be a progressively decreasing order, i.e., A > B > C > D.

[0106] Application Scenario 2 aims to generate simulated data showing significant differences in the activity of co-expression modules of the same gene. The specific steps are as follows: Four datasets with different cell numbers are generated using SERGIO simulation, corresponding to cell populations of 1000, 3000, 5000, and 10000 respectively. For each dataset, six cell populations (A, B, C, D, E, F) are simulated, with each cell population containing 100 genes, corresponding to four co-expression modules: three with similar activities (gene module 2, gene module 3, and gene module 4) and one with significantly different activities across the six cell populations (gene module 1). Assuming the drug's target gene is located in gene module 1, the actual drug response order across the six cell populations should be that A, B, and C rank significantly higher than D, E, and F, i.e., (A, B, C) > (D, E, F).

[0107] Application Scenario 3 aims to generate simulated data showing an ordered decrease in the activity of four co-expression modules of genes. The specific steps are as follows: Four datasets with different cell populations are generated using SERGIO simulation, corresponding to cell populations of 1000, 3000, 5000, and 10000 respectively. For each dataset, four cell populations (A, B, C, D) are simulated, each containing 100 genes, corresponding to four co-expression modules. In A, gene module 1 has significantly higher activity than the other three gene modules; in B, gene module 2 has significantly higher activity than the other three gene modules; in C, gene module 3 has significantly higher activity than the other three gene modules; and in D, gene module 4 has significantly higher activity than the other three gene modules. Assuming the drug's target genes are located in gene module 1 of A, gene module 2 of B, gene module 3 of C, and gene module 4 of D, the actual drug response order should be A, B, C, and D, respectively.

[0108] (2) Simulated data testing:

[0109] For disease samples driven by single or multiple pathogenic genes that interact with other genes, these genes often form a gene co-expression module. Therefore, three different types of simulated data and three application scenarios were designed to test the predictive performance of scRank.

[0110] The results of application scenario 1 show that scRank can accurately sort the four cell populations in the same order as the actual situation. Moreover, as the number of cells increases, the percentage of cell population A with the highest ranking is greater than the percentage of other cell populations.

[0111] The results of application scenario 2 show that scRank can distinguish between the "high group" and "low group" in gene module 1. The ranking difference is consistent with the actual situation. Moreover, as the number of cells increases, the overall ranking of A, B, and C is significantly higher than that of D, E, and F.

[0112] Results from application scenario 3 show that for perturbations of specific module genes, the scores of the scRank-predicted responding cell populations are higher than any other cell type, and under most conditions, the target cell population has an average scRank ranking of first. All these results demonstrate that scRank can accurately rank drug-responsive cell populations.

[0113] Case Study 1: Mouse Medulloblastoma.

[0114] Mouse medulloblastoma tumors and their scRNA-seq data after vismodegib treatment were obtained from the GEO platform, accession number GSE129730. The matrix data were standardized according to the authors' workflow (https: / / github.com / ben-babcock / Gershon_single-cell) and included 10 cell types: tumor cells (Node A, Node B, Node C, Node D), neurons, astrocytes, microglial cells, endothelial cells, fibroblasts, and oligodendrocytes. This study confirmed that vismodegib targets Smo, and that Node B showed sensitivity to vismodegib, while Node A showed resistance.

[0115] like Figure 4 As shown, Figure 4 This is a schematic diagram illustrating the application of scRank in mouse medulloblastoma data. Figure 4 (a) in the figure is a schematic diagram of the prediction results of scRank in mouse medulloblastoma, where the target gene Smo of vismodegib and the scRNA-seq data of mouse medulloblastoma in the disease state are the inputs of scRank; Figure 4 (b) in the figure is a schematic diagram of cell populations that are sensitive to and resistant to vismodegib, as confirmed by the literature, namely Node_B and Node_A, respectively; Figure 4 (c) in the diagram is a schematic diagram of the gene regulatory network of Node_A and Node_B in the disease state; Figure 4 In the diagram, (d) represents the GRNs in the two tumor subtypes and their differences before and after drug administration. The value in the diagram represents the weight of the edge in the GRN. Figure 4 (e) in the figure is a schematic diagram of the GO enrichment results of gene modules M1 and M2; Figure 4 (f) in the diagram is a schematic diagram of differentially expressed genes in Node_A and Node_B under disease conditions; Figure 4 A schematic diagram of the pathways that are significantly activated in (g)Node_B and Node_A; Figure 4 (h) in the figure represents the prognostic analysis of medulloblastoma patients based on the classification criteria of Node_A and Node_B characteristic genes.

[0116] scRank predicted Node_B to be the most responsive cell population to vemodilution, and Node_A to be the least responsive, a conclusion perfectly consistent with the actual results. The study found that Smo was a highly co-expressed gene in the GRN of Node_B during the disease state, while Smo was an isolated node in the GRN of Node_A. This partially explains why Node_B is sensitive to vemodilution, while Node_A is resistant. Cluster analysis of the gene co-expression matrices of Node_B and Node_A revealed that the difference in gene co-expression before and after drug administration was significantly higher in Node_B than in Node_A, and the gene co-expression module M1, containing Smo, was also significantly enriched with cell cycle-related signaling pathways. Further comparison of differentially expressed genes between Node_A and Node_B under disease states using Gene Set Enrichment Analysis (GSEA) revealed that Node_B significantly overexpressed genes such as Hey1, Mdk, Sox9, and Hes6, which were significantly positively correlated with the regulation of cell differentiation, while Node_A significantly overexpressed genes such as Ube2c, Top2a, Tpx2, and Ccna2, which were significantly positively correlated with cell division.

[0117] The above analyses all indicate that Node A is more malignant than Node B. Therefore, bulk sequencing data of medulloblastoma from 178 patients were collected. Based on the characteristic genes of Node A and Node B, these 178 patients were divided into two groups: the first group consisted of patients with high and low Node A expression, and the second group consisted of patients with high and low Node B expression. Combined with prognostic data, it was found that patients with high Node B expression had a better prognosis than those with low Node B expression, while patients with high Node A expression had a worse prognosis than those with low Node A expression. This further confirms that Node A is more malignant than Node B and validates the accuracy of scRank in predicting drug-sensitive and drug-resistant cell populations.

[0118] Case Study 2: Severe Depression.

[0119] The scRNA-seq data samples of the dorsolateral prefrontal cortex and its treatment with the selective serotonin reuptake inhibitor fluoxetine from patients with severe depression (MDD) were obtained from the GEO platform, accession number GSE144136. This dataset contains 25 cell types: excitatory neurons (Ex_1, Ex_2, Ex_3, Ex_4, Ex_5, Ex_6, Ex_7, Ex_8, Ex_9, Ex_10), inhibitory neurons (Inhib_1, Inhib_2, Inhib_3, Inhib_5, Inhib_6, Inhib_7, Inhib_8), endothelial cells, microglia, oligodendrocyte precursor cells (OPC_1, OPC_2), astrocytes (Astro_2, Astro_3), and oligodendrocytes (Oligos_1, Oligos_3). Literature confirms that fluoxetine is a commonly used drug for treating MDD. It mainly works by inhibiting SLC6A4 of excitatory neurons, thereby inhibiting serotonin transporter protein, blocking serotonin reuptake, increasing serotonin availability, and thus achieving the effect of treating MDD.

[0120] like Figure 5 As shown, Figure 5 This is a schematic diagram illustrating the application of scRank in dorsolateral prefrontal cortex data from patients with severe depression. Figure 5 (a) shows the scRNA-seq data of patients with severe depression and a schematic diagram of the target SLC6A4 of the therapeutic drug fluoxetine. Figure 5 (b) in the figure is a schematic diagram of the scRank ranking results of all cell populations on the severe depression patient dataset; Figure 5 (c) in the diagram is a schematic diagram of the gene regulatory network and cluster analysis results of Ex_9 and Inhib_5 in the disease state, which includes a co-expression module of 4 genes; Figure 5 (d) in the diagram is a schematic diagram of the pathway enrichment results of gene co-expression modules M1, M2, M3, and M4; Figure 5 (e) in the diagram is a schematic diagram of the activity distribution of the M2 gene co-expression module in all cell populations in the MDD; Figure 5 (f) in the figure is a schematic diagram of the differentially expressed genes of Ex_9 in normal and MDD; Figure 5 (g) in the diagram is a schematic diagram of the pathways that are significantly enriched in Ex_9, where the vertical axis represents the GSEA pathway score; Figure 5 (h) in the figure represents the spatial distribution of different cell populations in the mouse forebrain in MDD; Figure 5 (i) in the figure is a spatial distribution map of different cell populations in the human dorsolateral prefrontal cortex in MDD.

[0121] scRank predicted that excitatory neurons Ex_9 were the most responsive to fluoxetine, while Inhib_5 were the least responsive. To further explain this result, cluster analysis was performed on the gene co-expression matrices of Ex_9 and Inhib_5 under disease conditions. The results showed that the genes used to construct GRNs could be divided into four co-expression patterns, and the co-expression module M2, containing the fluoxetine target SLC6A4, was precisely the module with the highest activity in Ex_9. Enrichment results also confirmed that the co-expressed genes in M2 were significantly associated with neuronal synapse-related pathways, including Monoamine Transport, Serotonergic Synapse, and Synaptic Vesicle Pathway, which is consistent with the mechanism of action of fluoxetine. These results indicate that the gene co-expression module M2, containing SLC6A4, plays a crucial role in fluoxetine treatment of MDD. Further focusing on the M2 gene co-expression module, it was found that the M2 activity in Ex_9 was much higher than in other cell types, and the M2 activity in Ex_9 was also significantly higher than in the normal group. Meanwhile, the highly expressed genes in Ex_9 during the disease state were also significantly associated with synapses. These results indicate that the M2 gene co-expression module highly expressed in Ex_9 is highly correlated with the pathogenesis of MDD and the mechanism of action of fluoxetine in treating MDD.

[0122] Studies have shown that fluoxetine primarily acts on deep neurons. Therefore, cell populations in MDD were mapped onto spatial transcriptomic atlases of the mouse forebrain and human dorsolateral prefrontal cortex using the CellTrek tool (10X Visium datasets). The results showed that the Ex9 cells targeted by fluoxetine are indeed mainly distributed in deep neurons of the brain, with the majority of Ex9 cells concentrated in layers L5 and L6. These results further confirm that fluoxetine targets the Ex9 cell population to exert its therapeutic effect, inhibiting SLC6A4 of Ex9 cells and blocking serotonin reuptake, thereby achieving the therapeutic effect on MDD.

[0123] Five real-world scRNA-seq datasets with drug perturbation were collected from public databases through literature retrieval. Each dataset contains matrix data for both disease and post-drug administration states, clearly showing positive results for drug response and negative results for no response. Therefore, these five real-world datasets were used to benchmark scRank. Figure 6 As shown, its performance was compared with existing commonly used strategies and methods. Figure 6 A comparison graph of scRank with other methods; Figure 6 (a) in the figure is a schematic diagram comparing the performance of scRank, Augur and Wilcox using simulated data; Figure 6(b) in the figure is a schematic diagram of five sets of scRNA-seq datasets with real drug perturbation; Figure 6 (c) in the figure is a schematic diagram comparing the performance of scRank, Augur, Wilcox, MAST, bimod and DESeq2 using simulated data; Figure 6 (d) in the figure represents the box plot of the prediction results from different methods.

[0124] (1) Collection of real data.

[0125] Five real-world drug-perturbed scRNA-seq datasets were collected, including mouse medulloblastoma tumors, human colorectal cancer epithelial cells, mouse fibrotic bone marrow, miniaturized intestinal organoids, and mouse diabetic nephropathy datasets. The specific collection process is as follows:

[0126] Mouse medulloblastoma tumors and their scRNA-seq data after vismodegib treatment were obtained from the GEO platform, accession number GSE129730. The matrix data were standardized according to the authors' workflow (https: / / github.com / ben-babcock / Gershon_single-cell) and included 10 cell types: tumor cells (Node A, Node B, Node C, Node D), neurons, astrocytes, microglial cells, endothelial cells, fibroblasts, and oligodendrocytes.

[0127] Human colorectal cancer epithelial cells and their scRNA-seq data after leucovorin and 5-fluorouracil administration were obtained from the GEO platform, accession number GSE155953. A total of 11 different epithelial cell types were included: Epithelial cell 1, Epithelial cell 2, Epithelial cell 3, Epithelial cell 4, Epithelial cell 5, Epithelial cell 6, Epithelial cell 7, Epithelial cell 8, Epithelial cell 9, Epithelial cell 10, and Epithelial cell 11.

[0128] The scRNA-seq data samples of mouse fibrotic bone marrow and other quinomod-treated cells were obtained from the GEO platform, accession number GSE1 56644. The data included eight cell types: mesenchymal stromal cells (MSC_1, MSC_2, MSC_3, MSC_4), non-myelinating Schwanncell precursors (SCP), myelinating Schwanncell precursors (SCP), osteoblastic lineage cells (OLC), and endothelial cells.

[0129] Miniaturized intestinal organoids and their scRNA-seq data after KPT-330 treatment were obtained from the Single-CellPortal platform (https: / / singlecell.broadinstitute.org, SCP1547). This dataset contains 10 cell types: stem cells (stem cells I, II, III), enterocytes, early enterocytes, enteroendocrine cells, Paneth cells, early Paneth cells, early secretory cells, and goblet cells.

[0130] The scRNA-seq data samples of diabetic nephropathy in mice and after dapagliflozin treatment were obtained from the GEO platform, accession number GSE181382. The data included 10 cell types: proximal tubule cells, B cells, neutrophils, macrophages, collecting duct cells, endothelial cells, ascending loop of Henle cells, distal convoluted tubule cells, proliferating proximal tubule cells, and T cells.

[0131] (2) Real data test scRank.

[0132] The prediction performance of scRank was tested on five real datasets: mouse medulloblastoma tumors, human colorectal cancer epithelial cells, mouse fibrotic bone marrow, miniaturized intestinal organoids, and mouse diabetic nephropathy datasets. The specific results are as follows:

[0133] The first set of data consisted of mouse medulloblastoma tumors and their scRNA-seq data after vismodegib treatment, with vismodegib targeting Smo. This literature confirmed that Node_B was sensitive to vismodegib, while Node_A was resistant. The results showed that scRank predicted Node_B to be the most responsive cell population to vismodegib, and Node_A to be the least responsive, perfectly consistent with the actual results.

[0134] The second set of data consisted of human colorectal cancer epithelial cells and scRNA-seq data after administration of leucovorin and 5-fluorouracil, with leucovorin and 5-fluorouracil targeting TOP1 and TYMS, respectively. This study confirmed that epithelial cells 9, 10, and 11 were most sensitive to the combined administration of leucovorin and 5-fluorouracil, while epithelial cell 3 exhibited resistance. The results showed that scRank predicted epithelial cells 9, 10, and 11 to be in the top 50%, while epithelial cell 3 was the least responsive cell population, largely consistent with the actual results.

[0135] The third set of data consists of scRNA-seq data from mouse fibrotic bone marrow and other cells after taquimod administration, with taquimod targeting S100a8 and S100a9. This literature confirms that the target cell populations of taquimod are MSC1 and MSC2. The results show that scRank predicts that MSC1 and MSC2 cells will rank in the top 25% of drug responses, which is largely consistent with the actual results.

[0136] The fourth set of data consists of miniaturized intestinal organoids and their scRNA-seq data after KPT-330 treatment, with KPT-330 targeting Xpo1. This study confirms that KPT-330 targets stem cells II and III. The results show that scRank predicts stem cell III as having the top-ranking drug response, perfectly consistent with the actual results.

[0137] The fifth set of data consists of scRNA-seq data from diabetic nephropathy in mice and their response to dapagliflozin treatment, with dapagliflozin targeting Slc5a2. This study confirms that the target cell population for dapagliflozin is proximal tubules (PT). The results show that scRank predicted the top-ranking drug response in PT, perfectly consistent with the actual results.

[0138] (3) Compare scRank with other methods using real data.

[0139] Performance tests were conducted on both simulated and real data. Considering that other scRNA-seq data require two states, the "SERGIO" Python package and predefined GRN modules were used to automatically generate gene expression profile matrices for specific cell populations. The specific steps are as follows: Four datasets with different cell numbers were simulated using SERGIO, corresponding to cell populations of 1000, 3000, 5000, and 10000 respectively. For each dataset, four cell populations (A, B, C, D) were simulated, each containing 100 genes, corresponding to four gene co-expression modules. These included three gene co-expression modules with similar activities (gene modules M2, M3, and M4) and one gene co-expression module with gradually decreasing activity across the four cell populations (gene module M1), serving as scRNA-seq data for disease states. For each dataset, four corresponding cell populations (A', B', C', D') are simulated, each containing 100 genes and corresponding to four gene co-expression modules: three modules with similar activities (M2, M3, and M4) and one module (M1) with activity equal to M1 in A. Therefore, assuming the drug's target gene is located in module M1, the actual drug response order of the four cell populations after virtual knockout should be progressively decreasing, i.e., A>B>C>D.

[0140] The most common strategy currently used is to determine drug-responsive cell populations based on the number of differentially expressed genes before and after drug administration. Therefore, scRank was compared with the four most commonly used methods for screening differentially expressed genes: Wilcox, MAST, bimod, and DESeq2. Augur, Wilcox, MAST, bimod, and DESeq2 can generally accurately predict most drug-responsive cell populations, but scRank's prediction results are significantly better than the other methods. It is worth noting that other methods cannot accurately predict insensitive negative results. This demonstrates that scRank can predict not only drug-responsive cell populations but also drug-unresponsive cell populations.

[0141] Example 3

[0142] In order to implement the method corresponding to Embodiment 1 above and achieve the corresponding functions and technical effects, a drug response cell population ranking system based on single-cell transcriptome data is provided below.

[0143] A drug response cell population ranking system based on single-cell transcriptome data, comprising:

[0144] The data and information acquisition module is used to acquire scRNA-seq data and drug target information for disease states with known cell population types.

[0145] The GRN construction module is used to construct a target gene regulatory network (GRN) for each cell population based on the scRNA-seq data and the co-expression relationships of genes.

[0146] The tpGRN construction module is used to virtually knock out the target in the GRN according to the drug target information, and construct the target gene regulatory network tpGRN after virtual knockout of the target for each cell population.

[0147] A manifold alignment module is used to train a low-dimensional representation of each network node in the GRN and the tpGRN through manifold alignment; the network node is a gene node.

[0148] The Euclidean distance calculation module is used to calculate the distance of each network node in the GRN and the tpGRN based on the low-dimensional representation using Euclidean distance.

[0149] The drug perturbation score generation module is used to score the drug response of different cell populations based on the distance, taking into account the changing trends of the drug target, 2-hop node and the edge of the 2-hop node in the GRN and the tpGRN, and to generate a drug perturbation score.

[0150] The sorting module is used to sort the drug response of the cell population according to the drug perturbation score and determine the drug response cell population sorting result; the drug response cell population sorting result is used to characterize the degree of response of the cell population to drug perturbation.

[0151] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the systems disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the descriptions are relatively simple; relevant parts can be referred to the method section.

[0152] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for ranking drug response cell populations based on single-cell transcriptome data, characterized in that, The method comprises the following steps: Obtaining scRNA-seq data and drug target information of known cell population types in disease states; According to the scRNA-seq data, a target gene regulatory network GRN of each cell population is constructed by using the co-expression relationship of genes; According to the drug target information, the target points in the GRN are virtually knocked out to construct a target gene regulatory network tpGRN of each cell population after virtual knockout of the target points; The low-dimensional representation of each network node in the GRN and the tpGRN is trained by manifold alignment; The network node is a gene node; Based on the low-dimensional representation, the distance of each network node in the GRN and the tpGRN is calculated by Euclidean distance; Based on the distance, the drug response of different cell populations is scored by comprehensively considering the drug target, 2-hop node and the change trend of the edge of the 2-hop node in the GRN and the tpGRN, and a drug perturbation score is generated; According to the drug perturbation score, the drug response of the cell population is sorted, and a drug response cell population sorting result is determined; The drug response cell population sorting result is used to represent the response degree of the cell population after drug perturbation.

2. The method of ranking cell populations by drug response based on single cell transcriptome data according to claim 1, wherein, According to the scRNA-seq data, a target gene regulatory network GRN of each cell population is constructed by using the co-expression relationship of genes, specifically comprising: Constructing an scRNA-seq expression matrix according to the scRNA-seq data; Randomly selecting δc cells from each cell type in the scRNA-seq expression matrix, extracting the gene expression matrix corresponding to the selected δc cells, and repeating S times; wherein δ is the selection ratio, and c is the number of cells; Feature selection is performed on the genes in the gene expression matrix to construct a gene set of the GRN; the gene set includes the top 2000 high-variation genes in the scRNA-seq data, the genes corresponding to the transcription factors in the AnimalTFDB database, and the target genes of the drugs in the DGIdb database; The low-dimensional representation of each gene in the gene set is determined by principal component analysis, and the regression coefficient between genes is estimated as the weight of the edge in the GRN; Defining any gene expression matrix, using principal component analysis to regress the target gene on multiple potential covariates, and constructing the relationship between the target gene and the regulatory gene; According to the relationship between the target gene and the regulatory gene and the weight, a gene adjacency matrix is constructed; A tensor component analysis method is used to combine multiple gene adjacency matrices into a third-order tensor χ of S×T×R; wherein S is the number of GRNs, T is the number of target genes, and R is the number of regulatory genes; Extract the main feature mode shared by multiple GRNs, decompose and approximate the third-order tensor χ, and generate a new third-order tensor; Based on the new third-order tensor, each weight is divided by the maximum absolute value in the absolute value of all weights to normalize, and a final adjacency matrix of each cell population is constructed; the final adjacency matrix is a target gene regulatory network GRN.

3. The method of ordering cell populations by drug response based on single cell transcriptome data of claim 2, wherein, The relationship between the target gene and the regulatory gene is Y = X'β' + ε; wherein Y is an expression vector of a target gene; X' is an expression vector of a regulatory gene corresponding to a plurality of principal components ranked in the front; β' is a regression coefficient of a potential covariate after transformation; and ε is a random error.

4. The method of ordering cell populations by drug response based on single cell transcriptome data of claim 2, wherein, extracting a main feature pattern shared by a plurality of the GRNs, decomposing and approximating the third-order tensor χ to generate a new third-order tensor, specifically comprising: extracting a main feature pattern shared by a plurality of the GRNs, using the formula decomposing and approximating the third-order tensor χ; where s r is a vector of factors r in the GRN number dimension; is an outer product operation; t r is a vector of factors r in the target gene dimension; r r is a vector of factors r in the regulator gene dimension; selecting the first five factors to represent the key components of the third-order tensor χ, and generating a new third-order tensor through cumulative mean valueization.

5. The method of ordering cell populations by drug response based on single cell transcriptome data of claim 2, wherein, According to the drug target information, the target points in the GRN are virtually knocked out to construct a target gene regulatory network tpGRN after virtual knockout of the target points of each cell population, specifically comprising: extracting the row where the drug target gene is located from the final adjacency matrix, and setting the corresponding value of the row where the drug target gene is located to 0 to virtually construct an adjacency matrix after drug perturbation; the adjacency matrix after drug perturbation is the tpGRN.

6. The method of ranking cell populations by drug response based on single-cell transcriptome data according to claim 5, wherein, training the low-dimensional representation of each network node in the GRN and the tpGRN through manifold alignment, specifically comprising: combining the final adjacency matrix and the adjacency matrix after drug perturbation into a joint matrix; projecting the joint matrix to a shared low-dimensional latent space through manifold alignment to determine the low-dimensional representation of each gene node in the GRN and the tpGRN, and interpreting the Euclidean distance between each matched gene node as the perturbation effect of the drug on the target cell type.

7. The method of claim 6, wherein the method further comprises: determining a drug response cell population for each cell in the single-cell transcriptome dataset; and ranking the drug response cell populations based on the drug response cell population of each cell in the single-cell transcriptome dataset. The joint matrix is: wherein W is a joint matrix; A' is the final adjacency matrix; A PB is the adjacency matrix after drug perturbation; λ is an adjustment parameter; and I is an identity matrix describing the correspondence between the gene nodes in the GRN and the tpGRN.

8. The method of ranking cell populations by drug response based on single-cell transcriptome data of claim 1, wherein, The drug perturbation score score is: wherein, D Tar is the distance of the target gene in the GRN and the tpGRN; W out is the weight of the edge of the gene node leaving in the GRN; Deg Tar is the degree of the target gene in the GRN; N is the total number of genes regulated by the target gene; n is the gene serial number; D n is the distance of the downstream gene in the GRN and the tpGRN; W in is the weight of the edge of the gene node entering in the GRN; Deg n is the degree of the downstream gene in the GRN.

9. A drug response cell population ordering system based on single cell transcriptome data, characterized in that, comprising: a data and information acquisition module for acquiring scRNA-seq data and drug target information under the disease state of known cell population types; a GRN construction module for constructing a target gene regulatory network GRN of each cell population according to the co-expression relationship of genes using the scRNA-seq data; a tpGRN construction module for virtually knocking out the target points in the GRN according to the drug target information to construct a target gene regulatory network tpGRN after virtual knockout of the target points of each cell population; a manifold alignment module for training the low-dimensional representation of each network node in the GRN and the tpGRN through manifold alignment; the network node is a gene node; an Euclidean distance calculation module for calculating the distance of each network node in the GRN and the tpGRN based on the low-dimensional representation through Euclidean distance; a drug perturbation score generation module for scoring the drug response of different cell populations based on the distance, comprehensively considering the drug target, 2-hop node and the change trend of the edge of the 2-hop node in the GRN and the tpGRN to generate a drug perturbation score; a sorting module for sorting the drug response of cell populations according to the drug perturbation score to determine a drug response cell population sorting result; the drug response cell population sorting result is used to represent the response degree of the cell population after drug perturbation.

Images