Screening method of characteristic gene set and application of screening method in drug screening

Abstract

The invention discloses a characteristic gene set screening method and application thereof in drug screening, the method comprises the following steps: performing enrichment analysis on a function annotation gene set by using a gene expression data set, screening a significant gene set from the function annotation gene set by adopting a self-setting threshold value successive progressive analysis method, and screening the significant gene set from the function annotation gene set; and carrying out gene optimization to obtain feature gene sets of various grades, including self-setting a threshold value to carry out gene optimization on leading subsets with significant enrichment results. The characteristic gene set screening method is used for drug screening, can screen out target candidates with a gene expression regulating function from thousands of candidate disturbance data, is used for screening immune checkpoint inhibitor synergistic drugs, can effectively screen out compounds or compositions for enhancing ICB anti-cancer capacity, and has the advantages that the screening method is simple and convenient, and the screening efficiency is high. The existing medicine screening mode and frame are jumped out, and the existing technical problems are solved. Experimental results prove the effectiveness and accuracy of the method, target drugs can be accurately and massively screened out, and a powerful means is provided for discovery and development of new drugs.

Description

Technical Field

[0001] The present invention belongs to the field of biomedicine, and relates to a characteristic gene set and a drug screening method, and specifically to screening characteristic gene sets and drug screening based on the method of screening characteristic gene sets, and especially to screening of synergistic drugs of immune checkpoint inhibitors based on the characteristic gene set screening method and a characteristic gene set for screening synergistic drugs of immune checkpoint inhibitors. Background Art

[0002] Research on cancer immunotherapy has progressed rapidly. The 2018 Nobel Prize in Physiology or Medicine was awarded to Professor Tasuku Honjo and Professor James Allison, the discoverers of immune checkpoints PD-1 and CTLA-4. At present, a considerable number of tumors have been approved for treatment with immune checkpoint inhibitors (ICB). ICB immunotherapy is considered to be a powerful and effective clinical treatment that can even cure cancer patients. However, most patients with solid tumors do not respond to immune checkpoint inhibitor treatment. It is generally believed that the average response rate is only about 20%, and in some cancers it is even as low as 10% (highly invasive cancer). Researchers want to improve the effect of immune checkpoint inhibitor treatment through reasonable combination therapy, that is, to improve the response rate of immune checkpoint inhibitor treatment, or directly improve the efficacy of immune checkpoint inhibitors in treating tumors. Immunotherapy, especially immune checkpoint blockade (ICB) therapy (including anti-PD-1, anti-PD-L1 and anti-CTLA-4 therapy), aims to activate and mobilize the immune system to eliminate malignant cells. Although ICB therapy is considered a revolutionary treatment with the potential to cure different types of cancer, low objective efficacy rates remain a common problem for most patients with solid tumors, that is, a considerable proportion of patients do not respond to ICB therapy. Low response rates and treatment resistance remain current problems that need to be addressed, which hinder the effective use of ICB therapy for most patients. Therefore, it is crucial to explore combined therapies to improve the survival benefit of more patients. Currently, many studies are dedicated to improving the efficacy of ICB through combination therapy, that is, combining it with other treatments such as chemotherapy and targeted therapy. Some drugs have been reported to be used in combination with ICB therapy, which mainly enhance the efficacy of ICB therapy by regulating epigenetic states, targeting signaling pathways, increasing the number of cytotoxic T cells, and anti-angiogenesis. In these studies, most of the proposed drugs were discovered by chance, or were improved after further attempts based on existing mechanisms, or were modified based on existing drugs. However, due to the large-scale existence of existing drugs and drug targets, coupled with the heterogeneity of tumors, it is difficult to efficiently obtain and evaluate effective anti-tumor drugs combined with ICB. Currently, although the combination of ICB and chemotherapy has received attention, one of the main difficulties is the efficient prediction of compounds that can enhance the anti-tumor ability of ICB. Therefore, developing a more effective method to screen potential ICB combination drugs is an important issue facing the current industry. Summary of the invention

[0003] The task of the present invention is to provide a method for screening a characteristic gene set. Another task of the present invention is to provide a drug screening method based on the screening method of the characteristic gene set, including a screening method for synergistic drugs of immune checkpoint inhibitors based on the screening method of the characteristic gene set. Another task of the present invention is to provide a characteristic gene set for screening synergistic drugs of immune checkpoint inhibitors.

[0004] The technical solution for realizing the present invention is:

[0005] The method for screening a characteristic gene set provided by the present invention comprises the following steps G0, G1, G2 and G3:

[0006] Step G0: Gene expression datasets of humans and / or other animals are collected from public databases or experiments, and the collected gene expression datasets are grouped and the number of groups (N) is recorded. The principle of grouping is: each gene expression dataset is a group, or the gene expression datasets measured in the same batch are a group, or grouped by disease category, that is, the gene expression datasets of the same disease are a group; each group of gene expression datasets has at least two identical or similar phenotypic groups related to the research or application of the screening feature gene set;

[0007] Step G1: Collect gene sets with functional annotations;

[0008] Step G2: Use the gene expression data sets grouped by relevant phenotypes obtained in Step G0 to perform enrichment analysis on the gene sets collected in Step G1. The gene enrichment analysis is performed between the two identical or similar relevant phenotype groups; record the results of the enrichment analysis, including the normalized enrichment score NES, the significance test value (Q) of the enrichment result, and the leading-edge subset (L) of the significant enrichment result in the enrichment analysis result of each gene expression data set. The significance test value Q of the enrichment result is any one of p-value, adjust p-value, and FDR; and record the number (N1) of data sets in each gene expression data set group obtained in Step G0 where the enrichment results of each gene set with functional annotation collected in Step G1 are significantly positive and / or the number (N1') of data sets where the enrichment results are significantly negative; the method for judging a significantly positive enrichment result is: set a self-defined hypothesis test threshold (Q1); in the enrichment results where the gene set NES (enrichment score) > 0, when the significance test value Q of the enrichment result is less than the self-defined hypothesis test threshold Q1 set, that is, Q < Q1, it is a significantly positive enrichment result; preferably, the value of the self-defined hypothesis test threshold Q1 is less than 0.25; further preferably, the value of Q1 is 0.05, 0.01, 0.001, or 0.0001; the method for judging a significantly negative enrichment result is: set a self-defined hypothesis test threshold (Q1'); in the enrichment results where the gene set NES < 0, when the significance test value Q of the enrichment result is less than the self-defined hypothesis test threshold Q1', that is, Q < Q1', it is a significantly negative enrichment result; preferably, the value of the self-defined hypothesis test threshold Q1' is less than 0.25, and further preferably, the value of Q1' is 0.05, 0.01, 0.001, or 0.0001;

[0009] Step G3: Includes Step (A) and Step (B):

[0010] (A) Setting a screening threshold for a significant gene set in the enrichment result: including a screening threshold (X) for up-regulating a significant gene set and / or a screening threshold (X') for down-regulating a significant gene set; X and X' are independently valued, and the maximum value is equal to the number N of groups of the gene expression data set obtained in step G0, preferably X and / or X' are three quarters to one quarter of the value of the number N; further preferably X and / or X' are three quarters, one half, one third or one quarter of the value of the number N; X and X' are the same or Not the same; Screening of significantly up-regulated and / or down-regulated gene sets: The method for screening significantly up-regulated gene sets is: if the enrichment results of a gene set in more than or equal to X gene expression data sets are positively significant, that is, N1≥X, then the gene set is retained, and the retained gene set is a significantly up-regulated gene set; The method for screening significantly down-regulated gene sets is: if the enrichment results of a gene set in more than or equal to X' gene expression data sets are negatively significant, that is, N1'≥X', then the gene set is retained, and the retained gene set is a significantly down-regulated gene set;

[0011] (B) Optimizing the genes in each screened up-regulated and / or down-regulated significantly gene set: The method for optimizing the genes in each screened up-regulated significantly gene set is as follows: extracting the frontier subset L with significant enrichment results corresponding to each up-regulated significantly gene set in the enrichment analysis result, setting a gene optimization screening threshold (α), wherein the threshold α is at least 1 and at most equal to N1; detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each up-regulated significantly gene set according to the set gene optimization screening threshold α, retaining genes with an appearance greater than or equal to the threshold α, and discarding genes with an appearance less than the threshold α, to obtain an optimized gene set consisting of the retained genes, and the optimized gene set is obtained. The set is the up-regulated first characteristic gene set; the method for optimizing the genes in each screened out significantly down-regulated gene set is: extracting the frontier subset L with significant enrichment results corresponding to each significantly down-regulated gene set in the enrichment analysis result, setting a gene optimization screening threshold (α'), wherein the threshold α' is at least 1 and at most equal to N1'; according to the set gene optimization screening threshold α', detecting the appearance of each gene in the frontier subset L with significant enrichment results corresponding to each significantly down-regulated gene set, discarding genes with an appearance less than the threshold α', retaining genes with an appearance more than or equal to the threshold α', and obtaining an optimized gene set composed of the retained genes, and the optimized gene set is the down-regulated first characteristic gene set.

[0012] The at least two related phenotypic groups described in the above step G0 can be any one of the following ① to ③: ① "a group with good treatment effect after treatment" and "a group with no effect or little effect after treatment"; ② "a group of cancer patients receiving treatment" and "a group of normal people"; ③ "a group of cancer patients receiving treatment" and "a group of normal tissues or adjacent tissues of cancer patients receiving treatment". The at least two identical or similar phenotypic groups described in the above step G0 that are related to the purpose of the research or application of the screening feature gene set can specifically be: a "response" group and a "non-response" group after treatment with immune checkpoint inhibitors; the gene expression data set is a gene expression data set before and / or after treatment. The gene set with functional annotations described in the above step G1 is derived from any one or two or more of the following ① to ⑥: ① derived from the gene set with functional annotations officially provided by MSigDB; ② derived from the immune pathway gene set provided by the public resource library "ImmPort"; ③ derived from the annotated gene set provided in the KEGG public resource library; ④ derived from the annotated gene set provided in the WikiPathways public resource library; ⑤ derived from the gene set composed of gene markers corresponding to the same immune cell type, and the gene set composed of gene markers corresponding to the same immune cell type is specifically the marker gene set corresponding to the immune cells provided by the ImmuCellAI tool; ⑥ derived from the gene set composed of gene markers corresponding to immune cells with similar functions, and the gene set composed of gene markers corresponding to immune cells with similar functions is specifically the gene set composed of marker genes corresponding to killer cells provided by the ImmuCellAI tool.

[0013] The method for screening a characteristic gene set of the present invention further comprises a step G4 of screening a second characteristic gene set based on the above steps G1, G2 and G3.

[0014] Step G4: Screening the second feature gene set

[0015] First, the significant gene set to which the second characteristic gene set belongs is screened, including screening the up-regulated significant gene set and / or the down-regulated significant gene set: a screening threshold of the significant gene set to which the second characteristic gene set belongs is set, including a screening threshold (Y) for the up-regulated significant gene set and / or a screening threshold (Y') for the down-regulated significant gene set; Y and Y' are independently valued, Y is less than the screening threshold (X) for the up-regulated significant gene set to which the first characteristic gene set belongs, and Y' is less than the screening threshold (X') for the down-regulated significant gene set to which the first characteristic gene set belongs; the values ​​of Y and Y' are the same or different; the method for screening the up-regulated significant gene set is: according to the weights According to the enrichment analysis result of step G2 of claim 1, if the enrichment result of a certain gene set collected in step G1 is positively significant in more than or equal to Y and less than X data sets, then the gene set is retained, and the retained gene set is the significant gene set to which the second characteristic gene set is up-regulated; the method for screening the significantly down-regulated gene set is: according to the enrichment analysis result of step G2 of claim 1, if the enrichment result of a certain gene set collected in step G1 is negatively significant in more than or equal to Y' and less than X' gene expression data sets, then the gene set is retained, and the retained gene set is the significant gene set to which the second characteristic gene set is down-regulated;

[0016] Then, the genes in the significant gene set to which each up-regulated and / or down-regulated second characteristic gene set belongs are optimized:

[0017] The method for optimizing the genes in the significant gene set to which the screened up-regulated second characteristic gene set belongs is: extracting the frontier subset L with significant enrichment results corresponding to the significant gene sets to which each up-regulated second characteristic gene set belongs in the enrichment analysis result, setting a gene optimization screening threshold (β), wherein the threshold β is at least 1 and at most equal to Y, detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each up-regulated significant gene set according to the set gene optimization screening threshold β, retaining genes with more than or equal to the threshold β, and discarding genes with less than the threshold β, to obtain an optimized gene set composed of the retained genes, and the optimized gene set is the up-regulated second characteristic gene set;

[0018] The method for optimizing the genes in the significant gene set to which the screened down-regulated second characteristic gene set belongs is: extracting the frontier subset L with significant enrichment results corresponding to the significant gene sets to which each down-regulated second characteristic gene set belongs in the enrichment result; setting a gene optimization screening threshold (β'), wherein the threshold β' is at least 1 and at most equal to Y', detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each down-regulated significant gene set according to the set gene optimization screening threshold β', retaining genes with more than or equal to the threshold β', and discarding genes with less than the threshold β', thereby obtaining an optimized gene set composed of the retained genes, and the optimized gene set is the down-regulated second characteristic gene set.

[0019] The method for screening the characteristic gene set of the present invention further comprises a step G5 of screening a third characteristic gene set based on the above steps G1, G2, G3 and G4.

[0020] Step G5: Screening the third characteristic gene set

[0021] First, the significant gene set to which the third characteristic gene set belongs is screened, including screening the up-regulated significant gene set and / or screening the down-regulated significant gene set: setting the screening threshold of the significant gene set to which the third characteristic gene set belongs, including the screening threshold (Z) of the up-regulated significant gene set and / or the screening threshold (Z') of the down-regulated significant gene set; Z and Z' are independently valued, the threshold Z is less than the screening threshold Y of the up-regulated significant gene set to which the second characteristic gene set belongs, and Z' is less than the screening threshold Y' of the down-regulated significant gene set to which the second characteristic gene set belongs; the values ​​of Z and Z' are the same or different;

[0022] The method for screening the significantly up-regulated gene set is: according to the enrichment analysis result of step G2 of claim 1, if a gene set collected in step G1 has a positive and significant enrichment result in more than or equal to Z and less than Y data sets, then retain the gene set, and the retained gene set is the significant gene set to which the up-regulated third characteristic gene set belongs;

[0023] The method for screening down-regulated significantly gene sets is: according to the enrichment analysis result of step G2 of claim 1, if a gene set collected in step G1 has a negatively significant enrichment result in more than or equal to Z' and less than Y' gene expression data sets, then retain the gene set, and the retained gene set is the significant gene set to which the down-regulated third characteristic gene set belongs;

[0024] Then, the genes in the significant gene sets to which the up-regulated and / or down-regulated third characteristic gene sets are screened out are optimized:

[0025] The method for optimizing the genes in the significant gene sets to which each up-regulated third characteristic gene set belongs is as follows: extracting a frontier subset L with significant enrichment results corresponding to the significant gene sets to which each up-regulated third characteristic gene set belongs in the enrichment results, setting a gene optimization screening threshold (γ), wherein the threshold γ is at least 1 and at most equal to Z, detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each up-regulated significant gene set according to the set gene optimization screening threshold γ, retaining genes with an appearance greater than or equal to the threshold γ, and discarding genes with an appearance less than the threshold γ, to obtain an optimized gene set composed of the retained genes, wherein the optimized gene set is the up-regulated third characteristic gene set;

[0026] The method for optimizing the genes in the significant gene sets to which each down-regulated third characteristic gene set belongs is: extracting the frontier subset L with significant enrichment results corresponding to the significant gene sets to which each down-regulated third characteristic gene set belongs in the enrichment results, setting the gene optimization screening threshold (γ') to be retained, wherein the threshold γ' is at least 1 and at most equal to Z', detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each down-regulated significant gene set according to the set gene optimization screening threshold γ', retaining genes with more than or equal to the threshold γ', and discarding genes with less than the threshold γ', to obtain an optimized gene set composed of the retained genes, and the optimized gene set is the down-regulated third characteristic gene set.

[0027] The method for screening a characteristic gene set provided by the present invention can be used to screen drugs. Specifically, the drug screening method based on the screening method of the characteristic gene set provided by the present invention comprises the following steps:

[0028] Step 1: Obtain a corresponding characteristic gene set according to the method described in any one of claims 1 to 8. When multiple characteristic gene sets are obtained, merge the multiple characteristic gene sets into a parent characteristic gene set, which is a union of the multiple characteristic gene sets after deduplication;

[0029] Step 2: Collecting gene perturbation data of the screened compound or composition acting on cell lines, tissues and / or living experimental animals; the gene perturbation data specifically refers to the data of gene expression change values ​​when the screened compound or composition is acted on cell lines, tissues and / or living experimental animals, and then using gene expression detection technology, such as expression profile chip, RNA-seq or L1000, etc., to measure or calculate the gene expression change value relative to when the screened compound or composition is not added, that is, to obtain the gene perturbation data corresponding to each screened compound or composition;

[0030] Step 3: Using the enrichment analysis method to calculate the enrichment score NES and the enrichment significance test value P of the characteristic gene set or the parent characteristic gene set in each screened compound or composition perturbation data for the obtained gene perturbation data, wherein P is p-value, adjusted p-value or FDR;

[0031] When the screened target drug is an enhancer of the characteristic gene set or the parental characteristic gene set, the compound or composition corresponding to the screened enrichment analysis result meets the following conditions to be the target drug: NES>0 and P<P1, where P1 is a self-set threshold value of the enrichment significance test value. Preferably, P1 is less than 0.25. Further preferably, the value of P1 is 0.05, 0.01, 0.001 or 0.0001. When the screened target drug is an inhibitor of the parental characteristic gene set, the compound or composition corresponding to the screened enrichment analysis result meets the following conditions to be the target drug: NES<0 and P<P2, where P2 is a self-set threshold value of the enrichment significance test value. Preferably, P2 is less than 0.25. Further preferably, the value of P2 is 0.05, 0.01, 0.001 or 0.0001.

[0032] The screening method of the synergistic drug of the immune checkpoint inhibitor based on the characteristic gene set screening method provided by the present invention includes the following steps:

[0033] Step (1), screening the characteristic gene set, screening the characteristic gene set according to the method described in any one of claims 3-8 to obtain the corresponding characteristic gene set;

[0034] Step (2), collecting relevant gene expression data sets and grouping: collecting the gene expression data sets of the combined treatment of immune checkpoint inhibitors, and grouping the collected gene expression data sets of the combined treatment of immune checkpoint inhibitors into at least any two of the following four groups: <1> control group, <2> single use of synergistic treatment group, <3> single use of immune checkpoint inhibitor group, <4> combined use of immune checkpoint inhibitor and synergistic treatment group;

[0035] Performing enrichment analysis grouping on the gene expression data of each group of gene expression data sets to obtain one or two or more than two groups of the following 5 enrichment analysis groups: ① single use of synergistic treatment group vs. control group; ② single use of immune checkpoint inhibitor group vs. control group; ③ combined use of synergistic treatment group and immune checkpoint inhibitor group vs. single use of synergistic treatment group; ④ combined use of synergistic treatment group and immune checkpoint inhibitor group vs. control group; ⑤ combined use of immune checkpoint inhibitor treatment and synergistic treatment group vs. single use of immune checkpoint inhibitor group;

[0036] Recording the number of enrichment groups (H), that is, obtaining H enrichment groups, performing enrichment analysis on the obtained enrichment groups with each gene set in the characteristic gene set or the parental characteristic gene set to obtain the standardized enrichment score NES and the significance test value P0 and H enrichment results, where P0 is p-value, adjust p-value or FDR; recording the number (U) of those with significant enrichment and corresponding NES>0, and setting the NES values of these U enrichment results as positive sample data;

[0037] Step (III): Collect gene perturbation data of the screened compounds or compositions acting on cell lines, tissues, and / or living experimental animals. The gene perturbation data specifically refers to the data obtained by measuring or calculating the gene expression change values relative to the case without adding the screened compounds or compositions after the screened compounds or compositions act on cell lines, tissues, and / or living experimental animals, and then using gene expression detection techniques such as expression profiling chips, RNA-seq, or L1000, etc. That is, the gene perturbation data corresponding to each screened compound or composition is obtained; Combine the characteristic gene sets obtained in step (I) into a mother characteristic gene set, and use the enrichment analysis method to calculate the normalized enrichment score NES and the enrichment significance test value (P) of the mother characteristic gene set in the perturbation data of each screened compound or composition. The enrichment significance test value P is p-value, adjust p-value, or FDR;

[0038] When the screened target drug is an inhibitor of the mother characteristic gene set, the screened compounds or compositions corresponding to the enrichment analysis results meeting the following conditions are the target drugs: NES < 0 and P < P1, where P1 is a self-set enrichment significance test threshold; Preferably, the value of P1 is less than 0.25; Further preferably, the value of P1 is 0.05, 0.01, 0.001, or 0.0001; Set the NES value in the enrichment result corresponding to the perturbation data inhibited by the mother characteristic gene set as negative sample data, and use the positive sample data and negative sample data for machine learning modeling, and calculate the enrichment weight participated by each sub-gene set in the characteristic gene set. Preferably, machine learning methods such as LASSO regression or elastic net regression are used to obtain the perturbation score: R = X0 + a1X1 + a2X2 + a3X3 + ……; where R represents the perturbation score value, X0 represents the intercept in the machine learning calculation result, X1, X2, X3 …… represent the NES values in the GSEA results of the perturbation expression data of the characteristic gene set in the candidate compound, and a1, a2, a3 …… represent the weight coefficients corresponding to X1, X2, X3 ……; After obtaining the weight coefficients, score and rank each perturbation data in the perturbation data proposed as the target drug according to the perturbation score. The compound corresponding to the perturbation data with a higher score is more likely to be the target compound; Preferably, machine learning methods such as LASSO regression or elastic net regression are used for the machine learning modeling.

[0039] The synergistic therapy described herein refers to technical means that can be used in combination with immune checkpoint inhibitors to exert synergistic therapeutic effects, including but not limited to combined compounds, biological drugs, vaccines, radiotherapy, tumor-targeted drugs, gene knockout, gene overexpression and other technical means, that is, specifically, it can be combined compounds, biological drugs, vaccines, radiotherapy, tumor-targeted drugs, gene knockout and / or gene overexpression, etc.

[0040] The gene expression data set is RNA-seq raw data, expression profile chip raw data or gene expression matrix data obtained by upstream processing the raw data by common technical means in the field of bioinformatics, and the gene expression matrix data is Raw Counts, TPM, FPKM, RPKM gene expression matrix data or standardized gene expression matrix data. The public database can specifically be GEO, EBI, dbGAP and / or GSA.

[0041] The characteristic gene set for screening immune checkpoint inhibitor synergistic drugs provided by the present invention includes the first characteristic gene set C1-C4 and the second characteristic gene set M1-M2 obtained in the research of the present invention, and each characteristic gene set is specifically: the first characteristic gene set C1-C4:

[0042] C1: CD8A, HLA-DRA, HLA-DPA1, HLA-DQB1, CD74, IFNG, HLA-DRB5, HLA-DPB1, HLA-DQA1, HLA-DMA, HLA-DRB1, KLRC1, HLA-DQA2, CTSS, KL RD1, HLA-DOA, CIITA, CD8B, KIR2DL4, KLRC3, CD1B, HLA-DMB, CTSB, CD209, CD4, HLA-B, MICB, HLA-C, B2M, HLA-F, HLA-E, HLA-A, KLRC2;

[0043] C2: SH2D1A, IFNG, GZMB, LCP2, PTPN6, PIK3CD, ITGB2, LCK, FASLG, KLRD1, PRF1, HLA-E, CD247, CD48, CD 244, KIR2DL4, KLRC3, IFNB1, KLRC1, MICB, SH2D1B, FCER1G, HCST, HLA-B, HLA-C, HLA-A, TYROBP, KLRC2;

[0044] C3:LCK,CD3G,CD3D,ICOS,GRAP2,IFNG,LCP2,CD3E,PTPN6,PIK3CD,NFKBIA,CARD 11,CD247,ITK,PTPRC,PDCD1,CTLA4,PIK3R5,IL2,CD8B,CD8A,CD28,CD4,NFKBIE

[0045] C4:SIRPG,GPR171,CRTAM,GZMA,LAG3,CTSW,PRF1,NKG7,CCR5,C1QB,GZMB,G ZMH,LY9,CD7,LAX1,IL7R,ITK,IL2RB,LCP2,KLRG1,CELL,CD8B,CD8A,GNLY,

[0046] Removal range of snowmobiles M1-M2:

[0047] M1:CXCL11,APOBEC3A,CXCL9,IL2,ORM2,ISG20,TCHHL1,OASL,MX1,IL27,TNF,OAS1,CCL8,ISG15,BST2,MX2,CCL1, TLR7,MUC4,ORM1,DEFA3,CCR3,CHIT1,REG3G,C8G,CD40LG,CCR7,MPO,IDO1,IL22,CXCL10,CCL5,PAEP,GNLY,CD8A,C XCL13,CCL4,CCL7,CCL3,FASLG,FGR,APOBEC3H,HLA-B,MMP12,TLR8,APOBEC3G,IRF7,SLC29A3,HAMP,CD40,NOD2,C XCL2,CTSS,B2M,IRF5,IL10,MARCO,BPI,HCK,CYBB,IL6,CXCR1,DEFA4,DEFA1,S100A12,CAMP,PGLYRP1,CCL13,AZU1

[0048] M2:CD79A,CD19,CD79B,CARD11,CD72,INPP5D,PLCG2,PIK3CG,PTPN6,PRKCB,RAC2,VAV1,BTK,SYK,zPIK3CD,PIK3R5,NFKBIA,LYN,CR2,LILRB3,NFKBIE,JUN.

[0049] In specific use, some or all of the genes in the above gene sets C1, C2, C3, C4, M1 and M2 can be selected for screening immune checkpoint inhibitor synergistic drugs, that is, the characteristic gene set for screening immune checkpoint inhibitor synergistic drugs provided by the present invention can be a characteristic gene set composed of some or all of the genes selected from the above gene sets C1, C2, C3, C4, M1 and M2. The partial genes selected from the gene sets C1, C2, C3, C4, M1 and M2 can specifically be: at least 60% or 70% of the genes selected from any one of the gene sets C1, C2, C3, C4, M1 and M2, or at least 15 or 20 genes; or at least 40% or 50% of the genes selected from the union of any two or more gene sets in the gene sets C1, C2, C3, C4, M1 and M2, or at least 15 or 20 genes.

[0050] The immune checkpoint inhibitor described in the present invention can specifically be: anti-PD-1, anti-PD-L1 or anti-CTLA-4. The enrichment analysis described in the present invention is GSEA enrichment analysis; the enrichment result is the GSEA enrichment analysis result.

[0051] The present invention discloses a method for screening a characteristic gene set and its application in drug screening. The method uses a gene expression data set to perform enrichment analysis on a functional annotation gene set, adopts a self-set threshold and a progressive analysis method, screens out a significant gene set from the functional annotation gene set, and performs gene optimization on it, including setting a threshold to optimize the gene of the leading edge subset of the enrichment result significantly, and obtains characteristic gene sets of various levels. The characteristic gene set screening method is used for drug screening, including screening immune checkpoint inhibitor synergistic drugs, and can screen out target candidates with the function of regulating gene expression from tens of thousands of candidate compounds or composition perturbation data, and obtain substances that regulate gene expression and exert drug effects. The research experiment of the present invention proves the effectiveness and accuracy of the method, and can obtain a large number of effective candidate compounds at one time, jumping out of the existing drug discovery mode and framework, and can accurately and massively screen out target drugs, especially can efficiently obtain and evaluate drugs that can enhance the anti-tumor effect of ICB, and provides a powerful means for discovering and developing new drugs.

[0052] The screening method of the characteristic gene set provided by the present invention can screen the characteristic genes that cannot be screened by methods such as traditional differential genes. The annotation gene set can usually represent a certain biological process, and some genes of the annotation gene set have slight differences in expression between different phenotypic groups, and these differences have a unified law, but the traditional differential gene method cannot well identify these small but real changes in uniformity, because it must be omitted when screening according to the difference multiple or significance test value, and this omission leads to the inability to screen out these genes with the same slight changes in multiple similar data sets (all with groups of the phenotype of interest). And the method of the present invention can control the number of genes of the final characteristic gene set according to demand by setting multiple thresholds by itself, so that a gene set that meets the characteristics of multiple data sets can be screened.

[0053] The present invention uses bioinformatics methods to screen and predict ICB combination drugs. First, the core and secondary gene sets that respond to ICB therapy are determined by analyzing gene expression data. Then, based on the enrichment scores of the core and secondary gene sets, a new computational method, CM-Drug, is developed to predict the anti-tumor effect of ICB therapy combined with compounds. Through the CM-Drug method, we identified many potential ICB combination therapy compounds and verified the efficacy of the predicted drugs in melanoma and lung cancer mouse models. We also found a significant efficacy of a predictive drug, taltirelin, and explored its potential mechanism for combined use with ICB therapy. At present, most synergistic drug screening is based on targets or mechanisms, and it is difficult to jump out of the existing framework. There is currently no universal method for effectively screening immunotherapy combination drugs. In CM-Drug, a CM-Score is constructed through machine learning to predict and quantify the synergistic ability of compounds with ICB treatment. Existing human and mouse ICB combination data have confirmed the feasibility of our CM-Drug method. In addition, our mouse experiments on the predicted compounds also demonstrated the effectiveness of the CM-Drug method, which was confirmed in the repetition of different mouse experiments in different cancer types (such as melanoma and lung cancer) and datasets (LINCS Phase II and 2020 data). Our CM-Drug method fully utilizes the data of thousands of compounds provided in the LINCS L1000 database and breaks through the limitations of existing immunotherapy targets. In addition to predicting synergistic drugs in the LINCS database, our CM-Drug method can also be used to evaluate the therapeutic effects of ICB in combination with other compounds and even natural plant extracts. This will greatly accelerate the screening of ICB combination therapy drugs and greatly promote the clinical research of ICB therapy. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1It is a bubble chart visualization of the GSEA enrichment results (18 annotated gene sets) of 10 ICB clinical treatment data sets. The bubble with a check mark is a positively enriched result (NES>0). The size of the bubble represents the significance of the result. The value marked on the bubble is the normalized enrichment score NES (only the NES value corresponding to the adjusted p-value<0.05 is shown in the figure); see Example 1 for further explanation;

[0055] Figure 2 It is a schematic diagram for explaining the method used in step 3 of Example 1, taking Antigen Processing and Presentation (annotated gene set: antigen presentation) as an example;

[0056] Figure 3 The top 20 results of drug screening in the A375 cell line by the method described in Example 6, among which the three circled compounds are the compounds verified by in vivo mouse experiments in Example 7;

[0057] Figure 4 Schematic diagram of the administration method of Example 6, wherein Tumor bearing marked in the figure represents tumor bearing, that is, the 12th day before administration is the day when each group of mice bears tumor (-12 days), the arrows marked in the figure represent the days of administration, the day of administration is day 0, anti-PD-1 is administered 4 times, once every two days, Miglitol, MG-132 or Taltirelin are administered 6 times, once every one day, and ip in the figure represents intraperitoneal injection;

[0058] Figure 5 It is a tumor growth curve graph, specifically, the three compounds selected for the in vivo mouse test in Example 6 were combined with anti-PD-1, and the curve graph of the change in mouse tumor volume 9 days after treatment, each broken line represents the change in tumor volume growth of a corresponding mouse, and there are a total of 8 groups: 1, PBS (control group), 2, PBS + anti-PD-1 group, 3, Miglitol group, 4, Miglitol + anti-PD-1 group, 5, MG-132 group, 6, MG-132 + anti-PD1, 7, Taltirelin group, 8, Taltirelin + anti-PD-1 group

[0059] Figure 6This is a comparison of the tumor volume curves of mice in each group corresponding to the four groups of the efficacy experiment: PBS, anti-PD-1 antibody, Miglitol, and anti-PD-1+Miglitol (data are presented in the form of mean ± standard deviation). The asterisk represents that the difference in mouse tumor volume between the groups is extremely significant, that is, the efficacy of the combination of anti-PD-1 antibody and Miglitol is significantly better (* represents p-value < 0.05) than that of anti-PD-1 alone, which shows that Miglitol has the ability to enhance the treatment of anti-PD-1 antibody; in addition, Miglitol alone and anti-PD-1 alone are significantly different from the PBS group (control group);

[0060] Figure 7 This is a comparison of the tumor volume curves of mice in each group corresponding to the four groups of the efficacy experiment: PBS, anti-PD-1 antibody, MG-132, and anti-PD-1+MG-132 (data are presented in the form of mean ± standard deviation). The asterisk represents that the difference in mouse tumor volume between the groups is extremely significant, that is, the efficacy of anti-PD-1+MG-132 is significantly better (** represents p-value < 0.05) than that of anti-PD-1 antibody alone, indicating that MG-132 has a very significant ability to enhance the treatment of anti-PD-1 antibody. In addition, both MG-132 alone and anti-PD-1 alone have certain differences from the PBS group (control group), and the tumor volume is reduced relative to the PBS group, and the treatment effect of MG-132 alone after the 8th day is slightly better than that of the anti-PD-1 group alone, indicating that MG-132 alone has certain anti-tumor efficacy;

[0061] Figure 8 This is a comparison of the tumor volume curves of mice in each group corresponding to the four groups of the efficacy experiment: PBS, anti-PD-1 antibody, Taltirelin, and anti-PD-1+Taltirelin (data are presented in the form of mean ± standard deviation). The asterisk represents that the difference in mouse tumor volume between the groups is extremely significant, that is, the efficacy of the combination of anti-PD-1 antibody and Taltirelin is extremely significant (**** represents p-value < 0.0001) better than that of anti-PD-1 antibody or Taltirelin alone, indicating that Taltirelin has a very significant ability to enhance the treatment of anti-PD-1 antibody. In addition, Taltirelin alone and anti-PD-1 alone are significantly different from the PBS group (control group), and the treatment effect of Taltirelin alone after the 8th day is slightly better than that of the anti-PD-1 group alone, indicating that Taltirelin alone has a certain anti-tumor efficacy.

[0062] Fig. 9 This is a comparison of the survival curves of mice in the four groups of Taltirelin in the efficacy experiment of Example 7: PBS, anti-PD-1 antibody, Taltirelin, and anti-PD-1+Taltirelin. In the figure, "Log-rank (Mantel-Cox) test P value = 0.0004" represents the p-value of the Log-rank (Mantel-Cox) test = 0.0004, indicating that taltirelin significantly improves the survival of mice compared with the other groups;

[0063] Fig.10 : is a schematic diagram of the administration method of Example 8, wherein Tumor bearing marked in the figure represents tumor bearing, and the arrows marked in the figure represent the days of drug administration, that is, the day of administration is day 0, and the 7th day before administration is the day when each group of mice bears tumors (-7), anti-PD-1 is administered 4 times, once every two days, and taltirelin is administered 6 times, once every one day. On the 18th day, the tumors of each group of mice were cut off for subsequent experiments, and ip in the figure represents intraperitoneal injection;

[0064] Fig.11 It is a tumor growth curve diagram of the 4 groups of the pharmacodynamic experiment in Example 8 (the pharmacodynamic experiment of Taltirelin combined with anti-PD-1 was performed again), and each broken line represents the growth change of the tumor volume of a corresponding mouse; the 4 groups are: 1. PBS group (control), 2. anti-PD-1 single administration group, 3. Taltirelin single administration group, 4. anti-PD-1 + Taltirelin combined administration group;

[0065] Fig.12It is a comparison chart of the tumor volume curves of mice in each group corresponding to the four groups of the efficacy experiment in Example 8: PBS, anti-PD-1 antibody, Taltirelin, and anti-PD-1+Taltirelin (data are presented in the form of mean ± standard deviation) (data are presented in the form of mean ± standard deviation). The asterisk represents that the difference in mouse tumor volume between the groups is extremely significant, that is, the efficacy of the combination of anti-PD-1 antibody and Taltirelin is extremely significant (*** represents p-value <0.001, **** represents p-value <0.0001) better than anti-PD-1 antibody alone or Taltirelin alone, which shows that Taltirelin has a very significant ability to enhance anti-PD-1 antibody treatment. In addition, Taltirelin alone and anti-PD-1 alone are significantly different from the PBS group (control group), and the therapeutic effect of Taltirelin alone after the 8th day is slightly better than that of the anti-PD-1 group alone, which shows that Taltirelin alone has a certain anti-tumor efficacy;

[0066] Fig.13 In the pharmacodynamic experiment of Example 8, after the end of the pharmacodynamic experiment (18 days after tumor bearing), the tumors of all experimental mice in the four groups (PBS group, anti-PD-1 antibody group, Taltirelin group, anti-PD-1+Taltirelin group) were removed, and photos of all the removed tumors were taken;

[0067] Fig.14In the pharmacodynamic experiment of Example 8, at the end of the pharmacodynamic experiment, the tumors of all experimental mice in the four groups (PBS group, anti-PD-1 antibody group, Taltirelin group, anti-PD-1+Taltirelin group) were removed, and the tumor volume was measured and the relationship diagram of the tumor volume of each group of mice was drawn. The Kruskal-Wallis H method was used to test the difference in the data of the above four groups, and the difference among the four groups was significant. Pairwise tests were performed on the combination group and the other three groups: Mann-Whitney test was used for pairwise significance test, and the statistical differences between the two groups were significant (Adjust p-value<0.05). The False Discovery Rate BH method was used to correct them, and the correction results FDR were all less than 0.05, that is, the asterisks marked on the figure represent significant statistical differences between the two groups, indicating that after the end of the pharmacodynamic test, the final tumor volume of the anti-PD-1+Taltirelin group was significantly smaller than that of the single drug group (anti-PD-1 antibody group or Taltirelin group) and the solvent control group (PBS), which once again showed that Taltirelin combined with immune checkpoint inhibitor treatment has the effect of significantly reducing tumor volume;

[0068] Fig.15 This is a GSEA enrichment curve obtained by performing bulk RNA-seq sequencing on the excised tumor tissue after the efficacy experiment in Example 8, and then performing GSEA enrichment analysis on the sequencing results;

[0069] Fig.16 These are the top 20 compounds predicted in Example 10. The pills with shadow marks in the upper half of the figure are compounds that have been clearly reported to enhance the anti-tumor effect of ICB, and the compounds marked with circles are compounds that have no reports related to immunotherapy and are intended to be tested for efficacy in mice;

[0070] Fig.17 The dosing method and time of the mouse efficacy experiment in Example 11 are as follows: the day of the start of dosing is day 0, and each group of mice is tumor-bearing on day 9 before dosing. The mark in the bracket after the drug name in the figure represents the dosing method, where po represents oral gavage and ip represents intraperitoneal injection; anti-PD-1 is given once every three days (4 times of dosing, if the mouse dies during the period or the tumor volume exceeds the ethical limit, the dosing is stopped), and other drugs are given once every two days (6 times of dosing, if the mouse dies during the period or the tumor volume exceeds the ethical limit, the dosing is stopped);

[0071] Fig.18The results of the mouse efficacy experiment in Example 11 are shown in the figure. The figure shows the tumor inhibition curve. The asterisk represents the statistical test results (two-way ANOVA, followed by Bonferroni's multiple comparison test), * represents p-value < 0.05; ** represents p-value < 0.01; *** represents p-value < 0.001; **** represents p-value < 0.0001; the therapeutic effects of mice in the GSK429286A, CGP60474, camptothecin, wortmannin, AZD-7762 combined anti-PD-1 group are significantly better than those of mice in the anti-PD-1 group alone. The therapeutic effects of mice in the GSK429286A, CGP60474, camptothecin, AZD-7762 combined anti-PD-1 group are significantly better than those of the corresponding single drug treatment groups. DETAILED DESCRIPTION

[0072] Example 1: Screening of characteristic gene sets related to ICB treatment

[0073] 1. Gene expression data collection: The following table shows the clinical RNA-seq samples collected before ICB treatment, and the data comes from the NCBI GEO and dbGap databases; in order to include clinical treatment information on whether each patient responds to ICB treatment, these 10 data sets include 397 samples treated with anti-PD-1 alone, 91 samples treated with anti-PD-1 and anti-CTLA-4 at the same time, and 41 samples treated with anti-CTLA-4 alone. In order to avoid interference caused by different ICB treatments, we uniformly selected 400 samples treated with anti-PD-1 alone for further analysis. We defined complete response (CR) or partial response (PR) as the response (R) group, and stable disease (SD) or progressive disease (PD) as the non-response (NR) group. The relative phenotypes in this embodiment are: "response group after ICB treatment" and "non-response group after ICB treatment";

[0074]

[0075] 2. Collect annotated gene sets:

[0076] The initial gene set selection for screening: 17 immune pathways provided by the Immport website and cytotoxic cell marker genes provided by ImmuCellAI, a total of 18 gene sets (annotated gene sets).

[0077] 3. Obtain the significant difference markers between the "response" and "non-response" groups:

[0078] The total number of data sets collected is N=10, the conditional threshold X is set to 6, and the number of gene sets is 18; the hypothesis test value here uses the adjusted p-value, that is, Q is the adjusted p-value; the conditional threshold Q1 is set to 0.05, and the threshold α is set to 5.

[0079] As attached Figure 1 , in the enrichment results of all 18 gene sets, if the positive enrichment is significant, that is, the pathway that satisfies NES>0, p-value<0.05 and adjusted p-value<0.05 is marked with a "check mark" symbol. Since the number of data sets N=10, the conditional threshold X is set to 6, that is, if a gene set has 6 or more GSEA significant enrichments in 10 data sets, then the gene set is retained. As shown in the figure, the four data sets framed by solid lines: Antigen Processing and Presentation, TCR Signifying Pathway, Natural Killer Cell Cytotoxicity, Cytotoxicity of ImmuCellAI, these four gene sets are retained for further screening.

[0080] After retaining the four gene sets, the genes in these four gene sets were further screened. Taking Antigen Processing and Presentation as an example, the specific screening method is as follows Figure 2 As shown, Figure 2The bubble chart was converted into a more detailed visual curve chart. Antigen Processing and Presentation was significantly enriched in 7 data sets (a total of 10). The curves of the 7 data sets with significant enrichment were marked with a graphic "check mark", and the 3 data sets with insignificant enrichment were marked with a graphic "cross". The leading edge subsets of the 7 data sets with significant enrichment were extracted to obtain 7 leading edge subsets numbered 1-7. If a gene appears in 5 or more leading edge subsets, then the gene is retained. If it appears in less than 5 leading edge subsets, then the gene is discarded, and finally the optimized Antigen Processing and Presentation gene set is obtained. Through the above optimization, we screened out the gene set with the strongest signal, that is, the characteristic gene set screened by this method is the most significant gene set between the relevant phenotypes (the relevant phenotypes in this embodiment are "responding patients" and "non-responding patients" after receiving anti-PD-1 treatment, and the gene expression data are data before receiving treatment). Using these optimized gene sets to screen drugs greatly improves the accuracy of drug screening.

[0081] Similarly, the other three gene sets TCR Singnaling Pathway, Natural Killer Cell Cytotoxicity, and Cytotoxicity of ImmuCellAI were also screened in this way, and finally the gene sets C1-C4 were obtained as shown in the following table:

[0082]

[0083] The four gene sets shown in the above table are called the first characteristic gene sets, including four sub-gene sets C1-C4, which are the characteristic gene sets for target screening.

[0084] Example 2: Drug screening based on the characteristic gene set of Example 1

[0085] Collect perturbation data: Based on the perturbation data of a large number of small molecule compounds provided by the LINCS database, the characteristic gene set screened in Example 1 includes four characteristic gene sets C1-C4. The perturbation data (gene expression data change value) of all compounds or compositions are subjected to GSEA enrichment analysis on the four gene sets respectively, and the results of the four gene set enrichment results that simultaneously meet NES>0 and adjust p-value<0.05 are recorded as target candidate drugs.

[0086] Example 3: Further drug screening based on Example 2

[0087] A large number of candidate drugs were obtained in Example 2, including many drugs that have been reported to be proven to synergistically enhance the anti-tumor ability of ICB. However, since there are hundreds of target candidate drugs, further methods are needed to rank the hundreds of target candidate drugs. The characteristic gene set screened by Example 1 includes four characteristic gene sets C1-C4. Subsequently, the union of the four gene sets C1-C4 after deduplication is called the parent gene set, and all the perturbation data are used to perform GSEA enrichment analysis on the parent gene set. The results that meet NES>0 and adjust p-value<0.05 are the target candidate drugs. The results are sorted from large to small according to the NES value of GSEA, and the top ten drugs are as follows:

[0088]

[0089] Example 4: Adding a step of screening gene sets based on Example 1

[0090] On the basis of Example 1, the following steps are also included: In addition to the optimized four gene sets of TCR Singnaling Pathway, Natural Killer Cell Cytotoxicity, and Cytotoxicity of ImmuCellAI, we also screened out other gene sets. The more gene set conditions are defined, the more accurate the results will be in theory. The specific method is as follows:

[0091] Set the threshold Y=4 (according to Example 1, the number of data sets N=10, the conditional threshold X is set to 6, and the number of gene sets is 18; the hypothesis test value here uses the adjusted p-value, that is, Q is the adjusted p-value; the conditional threshold Q1 is set to 0.05, and the threshold β is set to 5)

[0092] As attached Figure 1 In the enrichment results of the 18 gene sets, if the positive enrichment is significant, that is, the pathways that meet NES>0, p-value<0.05 and adjust p-value<0.05 are marked with a "check mark" symbol. Since the number of data sets N1=10, the conditional threshold Y is set to 4, that is, if a gene set has significant GSEA enrichment in more than 4 and less than 6 of the 10 data sets, then the gene set is retained. Figure 1 As shown, the two data sets framed by dotted lines: Antimicrobials and BCR Signaling Pathway gene sets were retained for further screening.

[0093] After the two gene sets are retained, the genes in the two gene sets are further screened. Similar to the method in Example 1 ( Figure 2 ), Antimicrobials and BCR Signaling Pathway are significantly enriched in 4 data sets (a total of 10), and each significantly enriched data set will get a leading edge subset. If a gene appears in 3 or more leading edge subsets, then the gene is retained. If it appears in less than 3 leading edge subsets, then the gene is discarded. Finally, the optimized Antimicrobials and BCR Signaling Pathway gene sets are obtained, marked as M1 and M2 respectively. As shown in the following table

[0094]

[0095] The two gene sets shown in the above table are called the second characteristic gene sets, including two sub-gene sets M1-M2.

[0096] Example 5: Drug screening based on the characteristic gene sets screened in Example 1 and Example 4:

[0097] Collect perturbation data: The perturbation data of a large number of small molecule compounds provided by the LINCS database are used as the basis for drug screening. The LINCS Phase I and Phase II data contain a total of about 300,000 perturbation data of small molecule compounds. All perturbation data are used to perform GSEA enrichment analysis on 6 gene sets: optimized TCR Singnaling Pathway, Natural Killer Cell Cytotoxicity, Cytotoxicity of ImmuCellAI, Antimicrobials, and BCR Signaling Pathway. The results of these 6 gene sets that simultaneously meet NES>0, p-value<0.05, and adjust p-value<0.05 are recorded as target candidate drugs.

[0098] Example 6: Based on the method for screening drugs in Example 5, a screening step is added to further screen drugs

[0099] There are hundreds of candidate drugs in the screening results obtained in Example 5, and further steps are needed to rank them. The drugs with higher rankings are more likely to be the drugs that are most likely to cooperate with ICB to enhance the anti-tumor ability. The further screening steps of this example are as follows:

[0100] First, a gene expression data set of immune checkpoint inhibitor synergistic drug use is collected from a public data set. The gene expression data set of immune checkpoint inhibitor combination drug use includes one or more of the three groups of tumor control group and synergistic drug single use group, anti-PD-1 single use group, and anti-PD-1 combined with synergistic drug use group; the collected gene expression data can be: RNA-seq raw data, expression profile chip raw data, or gene expression matrix data obtained by upstream processing of the raw data by common technical means in the field of bioinformatics. The gene expression matrix data can be Raw Counts, TPM, FPKM, RPKM and other gene expression matrix data or standardized gene expression matrix data; the data accession numbers of the data sets collected from GEO for this modeling are: GSE130324, GSE149825, GSE126988, PRJNA477854, GSE146834, GSE148947, GSE149479, GSE137636, GSE132004, GSE79485. The above data sets all contain gene expression data for combined immunotherapy, and the drugs involved are as follows:

[0101]

[0102] The gene expression sequences of the immune checkpoint inhibitor synergistic drug data set collected above are calculated one by one for the 6 characteristic gene sets screened out in Example 4 (the 6 characteristic gene sets need to be converted into corresponding mouse gene sets by homologous genes) by the GSEA enrichment method to obtain the GSEA enrichment score ES, the standardized GSEA enrichment score NES, and the significance test result, i.e., p-value or adjusted p-value; each group of combination drug data includes: a tumor control group of the combination drug group and one or more of the three groups of the synergistic drug group alone, the anti-PD-1 group alone, and the combination of anti-PD-1 and synergistic drug. Then, the GSEA enrichment results of each group of combination drug data include at least one or more of the following four groups: synergistic drug group alone vs blank control group, anti-PD-1 group alone vs blank control group, combined synergistic drug group and anti-PD-1 group vs blank control group, combined synergistic drug group and anti-PD-1 group vs anti-PD-1 group alone; the GSEA enrichment results NES>0 and adjusted The NES value in the GSEA result with p-value<0.05 was used as the positive sample data for subsequent machine learning. The compound perturbation data (LINCS) was used to perform GSEA enrichment analysis on all perturbation data (gene expression change values), and the NES value of the GESA result that satisfied the adjusted p-value<0.05 and NES<0 was used as the negative sample data. The machine learning method was used to calculate the weight of each core gene set for the positive and negative samples, and the perturbation score was obtained: Y=0.4986+0.0969X1+0.0892X2+0.0307X3+0.0117X4+0.0124X5+0.0213X6, where Y represents the perturbation score value, 0.4986 represents the intercept in the machine learning calculation result, X1, X2, X3... represent the enrichment scores of the core gene set after GSEA standardization, and the coefficients of X1, X2, and X3 represent the weight of each gene set;

[0103]

[0104]

[0105] After obtaining the weight coefficient, each disturbance data in the target disturbance data is scored and sorted according to the disturbance score. The compound corresponding to the disturbance data with a higher score is more likely to be the target compound.

[0106] In order to verify whether the compounds screened by this method are really reliable, the results of the A375 cell line perturbation data are now displayed in the top 20 according to the above scores:

[0107]

[0108] Among the top 20 compounds with the highest scores, 6 have been reported to significantly cooperate with ICB in anti-tumor therapy, which indirectly proves the correctness of the screening method of the present invention, as follows:

[0109] Sequence 2: bortezomib (BRD-K50691590)

[0110] [1]Renrick AN, Thounaojam MC, Aquino M, et al. Bortezomib Sustains TCell Function by Inducing miR-155-Mediated Downregulation of SOCS1 and SHIP1[J]. Frontiers in Immunology, 2021,12:607044.

[0111] Mentioned in the article: These results support the approach that bortezomibcombined with other immunotherapies would lead to improved clinical outcomes by overcoming tumor-induced exhaustion or tolerance in CD8 + T cells. It is clearly reported that bortezomib can cooperate with other immunotherapies.

[0112] Ranking 4: Prochlorperazine (PCZ, prochlorperazine, BRD-K19352500)

[0113] Wrighton K H.Blocking endocytosis gives therapeutic antibodies aboost[J].Nature Reviews Drug Discovery, 2020,19(4).

[0114] The above-mentioned literature mentioned that PCZ enhanced the tumor growth inhibitory effect in BALB / c mice challenged with CT26.WT cells compared with ICB (avelumab) treatment alone.

[0115] Hui YC,Lima P,Cruz J,et al.Endocytosis Inhibition in Humans toImprove Responses to ADCC-Mediating Antibodies[J].Cell,2020,180(5):895-914.e27.

[0116] The highlights of this document describe that PCZ can synergistically enhance the efficiency of anti-tumor monoclonal antibodies.

[0117] Sort 12: AZD-7762

[0118] Peled M, Adam K, Mor A. Data on the identification of VRK2 as a mediator of PD-1 function. Data Brief. 2021May 21; 37:107168. doi: 10.1016 / j.dib.2021.107168. PMID: 34113705; PMCID: PMC8170101.

[0119] The above literature demonstrates that in a mouse model of MC38 tumor cells, combined treatment with AZD-7762 and anti-PD-1 enhanced the anti-tumor immune response compared with either treatment alone.

[0120] Ranking 16: LY-2157299 (Galunisertib, BRD-K78118466)

[0121] This drug has relevant reports in clinical trial research of the National Institutes of Health:

[0122] ClinicalTrials.gov Identifier:NCT02423343

[0123] A Study of Galunisertib(LY2157299) in Combination With Nivolumab inAdvanced Refractory Solid Tumors and in Recurrent or Refractory NSCLC,orHepatocellular Carcinoma

[0124] Study on the drug combined with Nivolumab in the treatment of advanced refractory solid tumors and recurrent or refractory NSCLC or liver cancer

[0125] Sorting 18: resiquimod (BRD-K76210423):

[0126] Nishii N, Tachinami H, Kondo Y, et al. Systemic administration of a TLR7 agonist attenuates regulatory T cells by dendritic cell modification and overcomes resistance to PD-L1 blockade therapy[J]. Oncotarget, 2018, 9(17): 13301-13312.

[0127] The above literature clearly reports: Our results suggest that systemic administration of low-dose resiquimod is useful as a companion drug to PD-1 / PD-L1 blockade therapy.

[0128] Sorting 20: propranolol (BRD-A10070317)

[0129] Kokolus K M, Zhang Y, Sivik J M, et al. Beta blocker use correlates with better overall survival in metastatic melanoma patients and improves the efficacy of immunotherapies in mice[J]. Oncoimmunology, 2017: 00-00.

[0130] The literature clearly reports: Treatment with propranolol alone provided a minimal, although statistically significant delay in tumor growth. However, propranolol administration significantly delayed tumor growth when combined with αPD-1 compared to either monotherapy. Propranolol with IL-2 also produced astatistically significant delay in tumor growth compared to IL-2alone, although the magnitude of improvement was less than observed in combination with αPD-1.

[0131] These reported drugs have confirmed the correctness of our method. It is worth noting that these reported drugs were obtained by chance or by experiments based on known pathways or targets, while the method provided by the present invention can jump out of the existing drug framework and relatively accurately find drugs that can enhance the anti-tumor effect of ICB from tens of thousands of compounds.

[0132] Example 7: Efficacy experiment to further prove that the drugs screened in Example 6 can enhance the anti-tumor effect of ICB

[0133] In addition to using a large number of reported drugs to prove the correctness of our screening method, we also want to discover new drugs through our method. We selected three compounds that ranked at the top of the results of Example 6 and have not been reported to have a combined immunotherapy enhancement effect, namely MG132, taltirelin, and miglitol, to conduct mouse experiments to verify whether they have a combined immunotherapy enhancement effect.

[0134] 1. Model construction and grouping

[0135] Source of reagents: Taltirelin, miglitol and MG-132 were purchased from Aladdin; anti-PD-1 was purchased from BioXCell

[0136] Experimental materials: C57BL / 6J male mice (6-8 weeks old), B16 / F10 subcutaneous tumor model

[0137] Grouping: PBS, Taltirelin, MG-132, Miglitol, PBS+anti-PD-1, Taltirelin+anti-PD-1, MG-132+anti-PD-1, Miglitol+anti-PD-1, 8 groups in total, 8 mice in each group. The administration method is as follows:

[0138]

[0139] PD-1 antibody administration: 100 μg per mouse per time, intraperitoneal injection, once every two days, for a total of four times.

[0140] 2. Pharmacodynamics experiment

[0141] After the start of drug administration, the tumor volume of mice was monitored every day. The values ​​in the table below are the tumor volume of mice (cubic millimeters). The values ​​in italics (asterisks) represent the death of mice midway, and only the data before death are available.

[0142] PBS group (control group)

[0143]

[0144]

[0145]

[0146] 3. Survival time

[0147] (1) Survival time data: The days in the table are calculated from the day the mice were injected with tumor cells, that is, the day of tumor bearing was day 0.

[0148]

[0149]

[0150] The following table shows the median survival time between the groups:

[0151]

[0152] Example 8: Pharmacodynamics experiment of Taltirelin combined with anti-PD-1

[0153] The results of Example 7 show that the three drugs screened in Example 6 all have the ability to synergize with anti-PD-1 in anti-tumor activities, among which Taltirelin has the best effect in enhancing anti-PD-1. A survival test was performed after the efficacy experiment in this example, and the efficacy test will be performed again in this example. After the efficacy test, the tumors of the mice were removed and weighed, and high-throughput sequencing (Bulk RNA-seq) was performed.

[0154] Experimental materials: C57BL / 6J male mice (6-8 weeks old); B16 / F10 subcutaneous tumor model (mouse melanoma cell line B16-F10 was purchased from the Cell Resource Center of Peking Union Medical College); Taltirelin was purchased from Aladdin Reagent Co., Ltd. with a specification of 10 mg; anti-PD-1 was purchased from BioXcell with a specification of 100 mg.

[0155] In this efficacy experiment, the following three groups, each with 8 mice, were used: PBS group, Taltirelin group (Taltirelin, a thyrotropin-releasing hormone analog), and PBS+anti-PD-1 group;

[0156] The following combination drug group was administered using 10 mice: Taltirelin + anti-PD-1.

[0157] The administration of Taltirelin and anti-PD-1 is as follows (see the dosing diagram for reference) Figure 5 ):

[0158]

[0159] 2. After the administration, the mouse tumors were removed and photographed. The tumor photos of each group of mice are attached to the instruction manual. Figure 7 The tumor weight of each group of mice was measured as follows (unit: grams). The statistical differences between the groups are shown in the appendix of the manual. Figure 8 .

[0160]

[0161]

[0162] 3. Statistical test

[0163]

[0164] Conclusion: The overall pharmacodynamics of the four groups was statistically significant. The statistical differences between the combination group and the other three groups were significant. The anti-tumor effect of Taltirelin combined with anti-PD-1 was significantly better than that of anti-PD-1 or Taltirelin alone. Figure 10-14 .

[0165] Example 9: Bulk RNA-seq sequencing

[0166] The mouse tumor tissue ( Fig.13 ) for RNA-seq sequencing, the mouse tumor tissue was removed and then quickly frozen in liquid nitrogen, then stored in a -80°C refrigerator, wrapped in a dry ice box and sent to Haplox for Bulk RNA-seq sequencing. Since the drug Taltirelin used in the experiment in Example 8 is a drug screened by the 6 characteristic gene sets (C1-C4, M1-M2) of Example 4, the sequencing results of Bulk RNA-seq were analyzed next. The results showed that the expression levels of the 6 gene sets in the combination drug (anti-PD-1+Taltirelin) group were extremely significantly improved compared with the two single drug groups (anti-PD-1 alone or Taltirelin alone) and the PBS control group; at the same time, the overall expression levels of the 6 gene sets in the responding samples were significantly higher than those in the non-responding samples (the tumor weight was less than the average tumor weight of the PBS group for the responding samples, and vice versa for the non-responding samples) ( Fig.15 ). This indicates that the combined drug therapy can significantly increase the overall expression level of the six characteristic gene sets, promote the immune response in the body, and thus better enable the body's own immune cells to play a role in killing tumor cells.

[0167] Example 10: New round of drug screening based on Example 4

[0168] This example uses the characteristic gene sets screened in Example 1 and Example 4 to screen drugs again. In this screening, we focused on the perturbation data of small molecule compounds in LINCS2020, and referred to the results of LINCS Phase I and LINCS Phase II. According to this method, a total of 6 characteristic gene sets (TCR singnaling Pathway, Natural Killer Cell Cytotoxicity, Cytotoxicity ofImmuCellAI, Antimicrobials, BCR Signaling Pathway) obtained by screening in Example 1 and Example 4 were subjected to GSEA enrichment analysis for all perturbation data, and the drugs corresponding to the perturbation data of these 6 gene sets that simultaneously meet NES>0, p-value<0.05, and adjust p-value<0.05 were found, which are the target candidate drugs. In order to verify the screening effect of this method in different cancers, the cell line for drug screening this time was set to the A549 human lung cancer cell line, and the mouse efficacy experiment will be used for verification in the subsequent Example 11. 536 candidate compounds were obtained in the A549 human lung cancer cell line in LINCS2020. In order to rank the candidate compounds, in this example, we calculated the total number of perturbation experiments of these 536 candidate compounds (536 compounds screened from the A549 human lung cancer cell line in LINCS2020) that met the screening criteria (6 feature gene sets NES>0, p-value<0.05, adjusted p-value<0.05) in the three phases of LINCS (Phase I, Phase II and 2020), focusing on the top 20 compounds (Tab.S5).

[0169]

[0170] Example 11: Verification of the efficacy of the compounds screened in Example 10 by mouse experiments

[0171] In order to determine the anti-tumor efficacy of the compounds screened in Example 10 in combination with anti-PD-1, we selected compounds that are easily accessible and have not been reported in immunotherapy-related studies for effect experiments. We further examined the top 20 candidate compounds ( Fig.16 ), 10 of which have been reported to have significant anti-tumor effects when combined with immunotherapy ( Fig.16), which directly illustrates the accuracy of the prediction of this method. Then we selected 6 of the remaining 10 compounds that have not been reported in immunotherapy for animal experiments: CGP60474, AZD-7762, camptothecin, GSK429286A, wortmannin and mycophenolic acid (MPA, BRD-K63750851).

[0172] In this study, the pharmacodynamics of the above six compounds were studied using C57BL / 6J mice bearing subcutaneous Lewis lung carcinoma (LLC) tumors.

[0173] Acquisition of experimental reagents: CGP60474, AZD-7762, camptothecin, GSK429286A, wortmannin and mycophenolic acid (MPA, BRD-K63750851) were all purchased from MedChemExpress.

[0174] Dosage and grouping: The drug was administered on the 9th day after the mice were tumor-bearing (at this time, the average tumor volume of each group of mice was about 100 mm 3 ). A total of 14 groups were divided into: 1. solvent control group 10% DMSO + 90% corn oil; 2. anti-PD-1 group alone; 3. GSK429286A group alone; 4. CGP60474 group alone; 5. camptothecin group alone; 6. wortmannin group alone; 7. AZD-7762 group alone; 8. mycophenolic acid group alone; 9. GSK429286A + anti-PD-1 group; 10. CGP60474 + anti-PD-1 group; 11. camptothecin + anti-PD-1 group; 12. wortmannin + anti-PD-1 group; 13. AZD-7762 + anti-PD-1 group alone; 14. mycophenolic acid + anti-PD-1 group.

[0175] There were 8 mice in each group. The dosage was calculated based on a 20g mouse. Specifically, the mice were dosed once every 2 days. The solvent control group was given 10% DMSO + 90% corn oil (Vehicle, i.e., the solvent control group, which was intraperitoneally injected with a dosing system of 100ul), GSK429286A (20mg / kg, oral gavage, dosing system of 200ul), CGP60474 (10mg / kg, intraperitoneally injected with a dosing system of 100ul), camptothecin (5mg / kg, intraperitoneally injected with a dosing system of 100ul), wortmannin (7mg / kg, intraperitoneally injected with a dosing system of 100ul), AZD-7762 (25mg / kg, intraperitoneally injected with a dosing system of 200ul), and mycophenolic acid (20mg / kg, intraperitoneally injected with a dosing system of 200ul). Acid (200 mg / kg, oral gavage, mice are calculated as 20 g each, and the dosing system is 200 ul), and a total of six injections; anti-PD-1 was injected intraperitoneally, and the dose of anti-PD-1 alone or in the combined group was one injection of anti-PD-1 every 3 days for each mouse, 100 μg each time. For specific dosing time and method, see Fig.17 Anti-PD-1 was administered 4 times in total. If the mouse died or the tumor volume exceeded the ethical limit during the period, the administration was stopped. Other drugs (compounds) were administered 6 times. If the mouse died or the tumor volume exceeded the ethical limit during the period, the administration was stopped.

[0176] Experimental results: mycophenolic acid (MPA, BRD-K63750851) did not enhance the anti-tumor ability of anti-PD-1, while the other five compounds GSK429286A, CGP60474, camptothecin, wortmannin, and AZD-7762 significantly improved the anti-tumor effect of anti-PD-1, as shown by a significant reduction in tumor volume in the combination group compared with the anti-PD-1 group alone ( Fig.18 ). These results demonstrate the effectiveness of the method of the present invention in different data and cancer types.

[0177] Those skilled in the art will readily appreciate that the above embodiments are not intended to limit the present invention, and any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for screening a characteristic gene set, comprising the following steps G0, G1, G2 and G3: Step G0: Collect gene expression data sets of humans and / or other animals from public databases or experiments, group the collected gene expression data sets and record the number of groups (N). The grouping principle is: each gene expression data set is a group, or gene expression data sets measured in the same batch are a group, or grouped by disease category, that is, gene expression data sets of the same disease are a group; each group of gene expression data sets has at least two identical or similar phenotype groupings related to the research or application for screening the characteristic gene set; Step G1: Collect gene sets with functional annotations; Step G2: Use the gene expression data sets of the relevant phenotype groupings obtained in Step G0 to perform enrichment analysis on the gene sets collected in Step G1. The gene enrichment analysis is performed between the two identical or similar relevant phenotype groupings; record the enrichment analysis results, including the normalized enrichment score NES, the significance test value (Q) of the enrichment result, and the leading-edge subset (L) of the significantly enriched result in the enrichment analysis result of each gene expression data set. The significance test value Q of the enrichment result is any one of p-value, adjust p-value, and FDR; and record the number (N1) of data sets with a positively significant enrichment result and / or the number (N1') of data sets with a negatively significant enrichment result of each gene set with functional annotations collected in Step G1 in the groupings of the gene expression data sets obtained in Step G0; The determination method for a positively significant enrichment result is: Set a self-set hypothesis test threshold (Q1); in the enrichment results where the gene set NES (enrichment score) > 0, when the significance test value Q of the enrichment result is less than the self-set hypothesis test threshold Q1 set, that is, Q < Q1, it is a positively significant enrichment result; The determination method for a negatively significant enrichment result is: Set a self-set hypothesis test threshold (Q1'); in the enrichment results where the gene set NES < 0, when the significance test value Q of the enrichment result is less than the self-set hypothesis test threshold Q1' set, that is, Q < Q1', it is a negatively significant enrichment result; Step G3: Includes Step (A) and Step (B) (A) Set the screening threshold for significantly enriched gene sets: Include the screening threshold (X) for up-regulated significantly enriched gene sets and / or the screening threshold (X') for down-regulated significantly enriched gene sets; X and X' each independently take values, and the maximum value is equal to the number of groups N of the gene expression data sets obtained in Step G0. Preferably, X and / or X' take values from three-quarters to one-quarter of the value of the number of groups N; further preferably, X and / or X' take values of three-quarters, one-half, one-third, or one-quarter of the value of the number of groups N; X and X' may take the same or different values; Screening for significantly up-regulated and / or down-regulated gene sets: The method for screening for significantly up-regulated gene sets is: if the enrichment results of a gene set in more than or equal to X gene expression data sets are positively significant, that is, N1≥X, then the gene set is retained, and the retained gene set is a significantly up-regulated gene set; the method for screening for significantly down-regulated gene sets is: if the enrichment results of a gene set in more than or equal to X' gene expression data sets are negatively significant, that is, N1'≥X', then the gene set is retained, and the retained gene set is a significantly down-regulated gene set; (B) Optimizing the genes in each screened up-regulated and / or down-regulated significantly gene set: The method for optimizing the genes in each screened significantly upregulated gene set is: extracting a frontier subset L with significant enrichment results corresponding to each significantly upregulated gene set in the enrichment analysis result, setting a gene optimization screening threshold (α), wherein the threshold α is at least 1 and at most equal to N1; According to the set gene optimization screening threshold α, the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each up-regulated significant gene set is detected, and genes with more than or equal to the threshold α are retained, and genes with less than the threshold α are discarded to obtain the optimized gene set composed of the retained genes, which is the up-regulated first characteristic gene set; The method for optimizing the genes in each significantly down-regulated gene set screened out is: extracting a frontier subset L with significant enrichment results corresponding to each significantly down-regulated gene set in the enrichment analysis result, setting a gene optimization screening threshold (α'), wherein the threshold α' is at least 1 and at most equal to N1'; According to the set gene optimization screening threshold α', the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each down-regulated significantly gene set is detected, and genes with less than the threshold α' are discarded, and genes with more than or equal to the threshold α' are retained to obtain the optimized gene set composed of the retained genes. The optimized gene set is the down-regulated first characteristic gene set.

2. The method according to claim 1, characterized in that In the method for judging whether the enrichment result is positively and / or negatively significant in step G2, the values ​​of the self-set hypothesis test thresholds Q1 and / or Q1' are: preferably, the values ​​of Q1 and / or Q1' are less than 0.25; further preferably, the values ​​of Q1 and / or Q1' are 0.05, 0.01, 0.001 or 0.0001; Q1 and Q1' are independently valued, and the values ​​of Q1 and Q1' are the same or different.

3. The method according to claim 1, characterized in that The at least two related phenotypic groupings described in step G0 are any one of the following ① to ③: ① "a group with good therapeutic effect after treatment" and "a group with no effect or little effect after treatment"; ② "a group of cancer patients receiving treatment" and "a group of normal people"; ③ "cancer tissue of cancer patients receiving treatment" and "normal tissue or adjacent tissue of cancer patients receiving treatment".

4. The method according to claim 1, characterized in that: The at least two identical or similar phenotypic groupings described in step G0 that are relevant to the purpose of research or application of the screening feature gene set are specifically: a "response" group and a "non-response" group after treatment with an immune checkpoint inhibitor; the gene expression data set is a gene expression data set before and / or after treatment.

5. The method according to claim 4, characterized in that The gene set with functional annotations described in step G1 is derived from any one or two or more of the following gene sets ① to ⑦: ① derived from the gene set with functional annotations officially provided by MSigDB; ② derived from the immune pathway gene set provided by the public resource library "ImmPort"; ③ derived from the annotated gene set provided by the KEGG public resource library; ④ derived from the annotated gene set provided by the WikiPathways public resource library; ⑤ derived from the annotated gene set provided by the Reactome public resource library; ⑥ derived from the gene set composed of gene markers corresponding to the same immune cell type; ⑦ derived from the gene set composed of gene markers corresponding to immune cells with similar functions.

6. The method according to claim 5, characterized in that The gene set composed of gene markers corresponding to the same immune cell type described in item ⑥ is the marker gene set corresponding to immune cells provided by the ImmuCellAI tool.

7. The method according to claim 5, characterized in that The gene set composed of gene markers corresponding to immune cells with similar functions as described in item 7 is specifically a gene set composed of marker genes corresponding to killer cells provided by the ImmuCellAI tool.

8. The method according to any one of claims 1 to 7, characterized in that: The following step G4 is also included: Step G4: Screening the second feature gene set First, the significant gene set to which the second characteristic gene set belongs is screened, including screening the up-regulated significant gene set and / or the down-regulated significant gene set: the screening threshold of the significant gene set to which the second characteristic gene set belongs is set, including the screening threshold (Y) of the up-regulated significant gene set and / or the screening threshold (Y') of the down-regulated significant gene set; Y and Y' are independently valued, Y is less than the screening threshold (X) of the up-regulated significant gene set to which the first characteristic gene set belongs, and Y' is less than the screening threshold (X') of the down-regulated significant gene set to which the first characteristic gene set belongs; the values ​​of Y and Y' are the same or different ; The method for screening the significantly up-regulated gene set is: if the enrichment result of a gene set collected in step G1 is positively significant in more than or equal to Y and less than X data sets, then the gene set is retained, and the retained gene set is the significant gene set to which the second characteristic gene set is up-regulated; The method for screening the significantly down-regulated gene set is: if the enrichment result of a gene set collected in step G1 is negatively significant in more than or equal to Y' and less than X' gene expression data sets, then the gene set is retained, and the retained gene set is the significant gene set to which the second characteristic gene set is down-regulated; Then optimize the genes in the significant gene set to which each up-regulated and / or down-regulated second characteristic gene set belongs: the method for optimizing the genes in the significant gene set to which the screened up-regulated second characteristic gene set belongs is: extract the frontier subset L with significant enrichment results corresponding to the significant gene set to which each up-regulated second characteristic gene set belongs in the enrichment analysis result, set a gene optimization screening threshold (β), the threshold β is at least 1 and at most equal to Y, detect the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each up-regulated significant gene set according to the set gene optimization screening threshold β, retain genes with more than or equal to the threshold β, discard genes with less than the threshold β, and obtain the optimized gene set composed of the retained genes, the The optimized gene set is the up-regulated second characteristic gene set; the method for optimizing the genes in the significant gene set to which the screened down-regulated second characteristic gene set belongs is: extracting the frontier subset L with significant enrichment results corresponding to the significant gene sets to which each down-regulated second characteristic gene set belongs in the enrichment result; setting a gene optimization screening threshold (β'), wherein the threshold β' is at least 1 and at most equal to Y', detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each down-regulated significant gene set according to the set gene optimization screening threshold β', retaining genes with more than or equal to the threshold β', and discarding genes with less than the threshold β', thereby obtaining an optimized gene set composed of the retained genes, and the optimized gene set is the down-regulated second characteristic gene set.

9. The method according to claim 8, characterized in that The following step G5 is also included: Step G5: Screening the third characteristic gene set First, the significant gene set to which the third characteristic gene set belongs is screened, including screening the up-regulated significant gene set and / or screening the down-regulated significant gene set: setting the screening threshold of the significant gene set to which the third characteristic gene set belongs, including the screening threshold (Z) of the up-regulated significant gene set and / or the screening threshold (Z') of the down-regulated significant gene set; Z and Z' are independently valued, the threshold Z is less than the screening threshold Y of the up-regulated significant gene set to which the second characteristic gene set belongs, and Z' is less than the screening threshold Y' of the down-regulated significant gene set to which the second characteristic gene set belongs; the values ​​of Z and Z' are the same or different; The method for screening the significantly up-regulated gene set is: if the enrichment result of a gene set collected in step G1 is positively significant in more than or equal to Z and less than Y data sets, then the gene set is retained, and the retained gene set is the significant gene set to which the up-regulated third characteristic gene set belongs; the method for screening the significantly down-regulated gene set is: if the enrichment result of a gene set collected in step G1 is negatively significant in more than or equal to Z' and less than Y' gene expression data sets, then the gene set is retained, and the retained gene set is the significant gene set to which the down-regulated third characteristic gene set belongs; Then, the genes in the significant gene sets to which the up-regulated and / or down-regulated third characteristic gene sets are screened out are optimized: The method for optimizing the genes in the significant gene sets to which each up-regulated third characteristic gene set belongs is as follows: extracting a frontier subset L with significant enrichment results corresponding to the significant gene sets to which each up-regulated third characteristic gene set belongs in the enrichment results, setting a gene optimization screening threshold (γ), wherein the threshold γ is at least 1 and at most equal to Z, detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each up-regulated significant gene set according to the set gene optimization screening threshold γ, retaining genes with an appearance greater than or equal to the threshold γ, and discarding genes with an appearance less than the threshold γ, to obtain an optimized gene set composed of the retained genes, wherein the optimized gene set is the up-regulated third characteristic gene set; The method for optimizing the genes in the significant gene sets to which each down-regulated third characteristic gene set belongs is: extracting the frontier subset L with significant enrichment results corresponding to the significant gene sets to which each down-regulated third characteristic gene set belongs in the enrichment results, setting the gene optimization screening threshold (γ') to be retained, wherein the threshold γ' is at least 1 and at most equal to Z', detecting the occurrence of each gene in the frontier subset L with significant enrichment results corresponding to each down-regulated significant gene set according to the set gene optimization screening threshold γ', retaining genes with more than or equal to the threshold γ', and discarding genes with less than the threshold γ', to obtain an optimized gene set composed of the retained genes, and the optimized gene set is the down-regulated third characteristic gene set.

10. Use of the method for screening characteristic gene sets according to any one of claims 1 to 9 in screening drugs.

11. A drug screening method comprising the following steps: Step 1: Obtain a corresponding characteristic gene set according to the method described in any one of claims 1 to 9. When multiple characteristic gene sets are obtained, merge the multiple characteristic gene sets into a parent characteristic gene set, which is a union of the multiple characteristic gene sets after deduplication; Step 2: Collecting gene perturbation data of the screened compound or composition acting on cell lines, tissues and / or living experimental animals; the gene perturbation data specifically refers to the data of the gene expression change value when the screened compound or composition is acted on cell lines, tissues and / or living experimental animals, and then using gene expression detection technology to measure or calculate the gene expression change value relative to when the screened compound or composition is not added, that is, obtaining the gene perturbation data corresponding to each screened compound or composition; Step 3: Using the enrichment analysis method to calculate the enrichment score NES and the enrichment significance test value P of the characteristic gene set or the parent characteristic gene set in each screened compound or composition perturbation data for the obtained gene perturbation data, wherein P is p-value, adjusted p-value or FDR; When the target drug to be screened is an enhancer of a characteristic gene set or a maternal characteristic gene set, the compound or composition corresponding to the screened compound or composition is a target drug if the enrichment analysis result satisfies the following conditions: NES > 0 and P < P1, where P1 is a self-set threshold value for the enrichment significance test; preferably, P1 is less than 0.25; more preferably, the value of P1 is 0.05, 0.01, 0.001 or 0.0001; when the target drug to be screened is an inhibitor of the maternal characteristic gene set, the compound or composition corresponding to the screened compound or composition is a target drug if the enrichment analysis result satisfies the following conditions: NES < 0 and P < P2, where P2 is a self-set threshold value for the enrichment significance test; preferably, P2 is less than 0.25; more preferably, the value of P2 is 0.05, 0.01, 0.001 or 0.0001.

12. A screening method for an immune checkpoint inhibitor synergistic drug, comprising the following steps: Step (1), screening a characteristic gene set, screening a characteristic gene set according to the method described in any one of claims 3-9 to obtain a corresponding characteristic gene set; Step (2), collecting relevant gene expression data sets and grouping: collecting gene expression data sets of immune checkpoint inhibitor combination therapy, and grouping the collected gene expression data sets of immune checkpoint inhibitor combination therapy into at least any two of the following four groups: <1> control group, <2> single-agent synergistic treatment group, <3> single-agent immune checkpoint inhibitor group, <4> combined immune checkpoint inhibitor and synergistic treatment group; Performing enrichment analysis grouping on the gene expression data of each group of gene expression data sets to obtain one group or two groups or more than two groups of the following 5 enrichment analysis groupings: ① single-agent synergistic treatment group vs. control group; ② single-agent immune checkpoint inhibitor group vs. control group; ③ combined synergistic treatment group and immune checkpoint inhibitor group vs. single-agent synergistic treatment group; ④ combined synergistic treatment group and immune checkpoint inhibitor group vs. control group; ⑤ combined immune checkpoint inhibitor treatment and synergistic treatment group vs. single-agent immune checkpoint inhibitor group; Recording the number of enrichment groupings (H), that is, obtaining H enrichment groupings, performing enrichment analysis on the obtained enrichment groupings with each gene set in the characteristic gene set or the maternal characteristic gene set to obtain a standardized enrichment score NES, a significance test value P0, and H enrichment results, where P0 is p-value, adjust p-value or FDR; recording the number (U) of significantly enriched ones with corresponding NES > 0, and setting the NES values of these U enrichment results as positive sample data; Step (III): Collect gene perturbation data of the screened compounds or compositions acting on cell lines, tissues, and / or live experimental animals. The gene perturbation data specifically refers to the data obtained by measuring or calculating the gene expression change values when the screened compounds or compositions act on cell lines, tissues, and / or live experimental animals, and then using gene expression detection techniques, that is, obtaining the gene perturbation data corresponding to each screened compound or composition; Combine the characteristic gene sets obtained in step (I) into a mother characteristic gene set, and use the enrichment analysis method to calculate the normalized enrichment score NES and the enrichment significance test value (P) of the mother characteristic gene set in the perturbation data of each screened compound or composition. The enrichment significance test value P is p-value, adjust p-value, or FDR; When the screened target drug is an inhibitor of the mother characteristic gene set, the screened compounds or compositions corresponding to the enrichment analysis results that meet the following conditions are the target drugs: NES < 0 and P < P1, where P1 is a self-set enrichment significance test threshold; preferably, the value of P1 is less than 0.25; more preferably, the value of P1 is 0.05, 0.01, 0.001, or 0.0001; set the NES value in the enrichment result corresponding to the perturbation data inhibited by the mother characteristic gene set as negative sample data, and use the positive sample data and negative sample data to perform machine learning modeling, and calculate the enrichment weight participated by each sub-gene set in the characteristic gene set. Preferably, use machine learning methods such as LASSO regression or elastic net regression to obtain the perturbation score: R = X0 + a1X1 + a2X2 + a3X3 +...; where R represents the perturbation score value, X0 represents the intercept in the machine learning calculation result, X1, X2, X3... represent the NES values in the GSEA results of the perturbation expression data of the characteristic gene set in the candidate compound, and a1, a2, a3... represent the weight coefficients corresponding to X1, X2, X3...; after obtaining the weight coefficients, score and sort each perturbation data in the perturbation data proposed as the target drug according to the perturbation score. The higher the score of the perturbation data, the more likely the corresponding compound is the target compound; preferably, use machine learning methods such as LASSO regression or elastic net regression to perform the machine learning modeling.

13. A characteristic gene set, characterized in that: The characteristic gene set is used to screen immune checkpoint inhibitor synergistic drugs, and the characteristic gene set includes some or all of the genes selected from the following gene sets C1, C2, C3, C4, M1, and M2; the gene sets C1, C2, C3, C4, M1, and M2 are respectively: C1:CD8A,HLA-DRA,HLA-DPA1,HLA-DQB1,CD74,IFNG,HLA-DRB5,HL A-DPB1,HLA-DQA1,HLA-DMA,HLA-DRB1,KLRC1,HLA-DQA2,CTSS,KL RD1,HLA-DOA,CIITA,CD8B,KIR2DL4,KLRC3,CD1B,HLA-DMB,CTSB, CD209,CD4,THE-B,MICB,THE-C,B2M,F-F,E-E,A,KLRC2 C2:SH2D1A,IFNG,GZMB,LCP2,PTPN6,PIK3CD,ITGB2,LCK,FASLG,KLRD1,PRF1,HLA-E,CD247,CD48,CD 244,KIR2DL4,KLRC3,IFNB1,KLRC1,MICB,SH2D1B,FCER1G,HCST,HLA-B,HLA-C,HLA-A,TYROBP,KLRC2 C3:LCK,CD3G,CD3D,ICOS,GRAP2,IFNG,LCP2,CD3E,PTPN6,PIK3CD,NFKBIA,CARD 11,CD247,ITK,PTPRC,PDCD1,CTLA4,PIK3R5,IL2,CD8B,CD8A,CD28,CD4,NFKBIE C4:SIRPG,GPR171,CRTAM,GZMA,LAG3,CTSW,PRF1,NKG7,CCR5,C1QB,GZMB,G ZMH,LY9,CD7,LAX1,IL7R,ITK,IL2RB,LCP2,KLRG1,CELL,CD8B,CD8A,GNLY M1: CXCL11, APOBEC3A, CXCL9, IL2, ORM2, ISG20, TCHHL1, OASL, MX1, IL27, TNF, OAS1, CCL8, ISG15, BST2, MX2, CCL1, TLR7, MUC4, ORM1, DEFA3, CCR3, CHIT1, REG3G, C8G, CD40LG, CCR7, MPO, IDO1, IL22, CXCL10, CCL5, PAEP, GNLY, CD8A, C XCL13, CCL4, CCL7, CCL3, FASLG, FGR, APOBEC3H, HLA-B, MMP12, TLR8, APOBEC3G, IRF7, SLC29A3, HAMP, CD40, NOD2, CX CL2, CTSS, B2M, IRF5, IL10, MARCO, BPI, HCK, CYBB, IL6, CXCR1, DEFA4, DEFA1, S100A12, CAMP, PGLYRP1, CCL13, AZU1; M2: CD79A, CD19, CD79B, CARD11, CD72, INPP5D, PLCG2, PIK3CG, PTPN6, PRKCB, RAC2, VAV1, BTK, SYK, zPIK3CD, PIK3R5, NFKBIA, LYN, CR2, LILRB3, NFKBIE, JUN.

14. The characteristic gene set according to claim 13, characterized in that: The partial genes selected from gene sets C1, C2, C3, C4, M1 and M2 are specifically: comprising at least 60% or 70% of the genes, or at least 15 or 20 genes, selected from any one of the gene sets C1, C2, C3, C4, M1 and M2; or It includes at least 40% or 50% of the genes, or at least 15 or 20 genes, selected from the union of any two or more gene sets among the gene sets C1, C2, C3, C4, M1 and M2.

15. The method according to claim 4, 12 or 13, characterized in that The immune checkpoint inhibitors mentioned are: anti-PD-1, anti-PD-L1 or anti-CTLA-4.

Images