Establishment and verification method of melanoma anti-PD-1 immunotherapy response predictive gene pair

Abstract

The invention discloses a method for establishing and verifying a melanoma anti-PD-1 immunotherapy response predictive gene pair. The method comprises the following steps: selecting a melanoma anti-PD-1 immunotherapy data set as a training set, forming a gene pair matrix by adopting a method based on a relative expression rank relationship of genes in a sample, and screening candidate predictive gene pairs; and on the basis of the candidate prediction gene pairs, a gene pair combination with the maximum sample coverage rate is searched by adopting a greedy algorithm, and the gene pair combination with the best prediction efficiency is screened as a prediction gene pair. Aiming at a to-be-predicted sample, calculating and predicting geneagt in a gene pair combination; and the gene logarithm of geneb. According to the voting rule, if the voting is geneagt; if the number of the gene pairs of the geneb is more than half, the to-be-predicted sample can be classified as a response sample, otherwise, the to-be-predicted sample is a non-response sample. And verification is carried out in sample verification sets of other different platforms.

Description

Technical Field

[0001] The present invention belongs to the technical field of melanoma anti-PD-1 immunotherapy drug response prediction, and specifically relates to a method for establishing and verifying a melanoma anti-PD-1 immunotherapy response prediction gene pair. Background Art

[0002] Melanoma is a highly invasive malignant tumor that often occurs in the skin, mouth, intestines, eyes and other parts of the body. It is easy to metastasize and difficult to treat. In recent years, anti-PD-1 immunotherapy has been used to treat patients with advanced melanoma and has achieved some success. However, due to the heterogeneity of tumors and individual patients, anti-PD-1 immunotherapy is not suitable for all patients. Some patients cannot benefit from this therapy and even experience accelerated cancer progression. At the same time, anti-PD-1 immunotherapy is expensive and has huge medical costs.

[0003] Studies have shown that the effect of anti-PD-1 immunotherapy is related to a variety of molecular markers (such as PD-L1, LAG-3, CXCR3, etc.), tumor mutation load, tumor lymphocyte infiltration level, new antigens and other factors. However, the actual application of these biomarkers is easily affected by factors such as expression level and batch effect, which makes cross-dataset and cross-platform verification difficult. These reasons are likely to affect the actual clinical application of biomarkers. Summary of the invention

[0004] In view of the fact that the effect of anti-PD-1 immunotherapy in the prior art is related to multiple molecular markers (such as PD-L1, LAG-3, CXCR3, etc.), tumor mutation load, tumor lymphocyte infiltration level, new antigens and other factors, but the actual application of these biomarkers is easily affected by factors such as expression level and batch effect, resulting in greater difficulty in cross-dataset and cross-platform verification, the present invention provides a method for establishing and verifying a gene pair for predicting the response to anti-PD-1 immunotherapy of melanoma, and provides a method for finding stable predictive drug response markers with clinical translation value in melanoma anti-PD-1 immunotherapy samples.

[0005] In order to solve the above technical problems, the present invention adopts the following technical solutions: A method for establishing a gene pair predicting response to anti-PD-1 immunotherapy for melanoma comprises the following steps: Step A: obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from a comprehensive gene expression database and a cancer genome atlas, and preprocessing the expression profile data and clinical information; The melanoma expression profile datasets from different detection platforms were divided into two data subsets, one melanoma anti-PD-1 immunotherapy dataset was used as a training set, and the other data subset was used as a validation set; Step B: Combine genes in the training set in pairs, compare the differences in relative expression rank relationships between gene pairs in progressive disease (PD) samples and partial response and complete response (PRCR) samples, screen gene pairs with significantly reversed relative expression rank relationships between genes, and identify candidate predictive gene pairs with potential ability to predict the response to anti-PD-1 immunotherapy for melanoma; Step C: Use a greedy algorithm to find the gene pair combination with the maximum sample coverage; Starting from each gene pair in the candidate prediction gene pair, use it as a seed and add the remaining gene pairs to the current combination in turn; after each addition, calculate the sample coverage of the new combination; this process continues until the sample coverage can no longer be increased by adding candidate gene pairs; the final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; Calculate the frequency of each characteristic gene pair in all combinations. The frequency reflects the importance of the characteristic pair. Select the best prediction gene pair as the final prediction gene pair. For a sample to be predicted, calculate the gene in the predicted gene pair combination a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR); the prediction gene pair is then verified using the validation set, and the obtained prediction gene pair is used as a marker for predicting the response to anti-PD-1 immunotherapy of melanoma. Further, the step A of obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from the gene expression comprehensive database and the cancer genome atlas, and preprocessing the expression profile data and clinical information specifically includes the following steps: Step A1: Download expression profile data from the GEO database and use the platform annotation file to map the probe ID to the gene ID; Step A2: Studies have shown that in the process of anti-PD-1 immunotherapy for melanoma, the subsequent efficacy of samples in the stable disease (SD) state is not stable, and they may develop into PD samples or PRCR samples in subsequent treatments. Therefore, SD samples are removed from the training set; Step A3: Remove low-expression genes, which are genes whose expression is missing or zero in more than 70% of sample genes.

[0006] Furthermore, in step A1, the expression profile data is downloaded from the GEO database, and the probe ID is mapped to the gene ID using the platform annotation file. The melanoma anti-PD-1 immunotherapy expression profile dataset GSE91061 (including samples before and during treatment) from the GPL9052 detection platform is used as a training set (GSE91061_Pre) and a validation set (GSE91061_On), respectively, and the other three sets of melanoma expression profile datasets are used as independent validation sets; the three sets of validation sets are melanoma expression profile data from TCGA, GSE115821, and GSE22153 of the detection platforms GPL11154, GPL18573, and GPL6102.

[0007] Furthermore, the genes in the training set are combined in pairs as described in step B, the differences in the relative expression rank relationships of the gene pairs in the PD samples and the PRCR samples are compared, the gene pairs with significantly reversed relative expression rank relationships between the genes are screened, and the candidate predictive gene pairs with potential ability to predict the response to anti-PD-1 immunotherapy of melanoma are determined, which specifically includes the following steps: Step B1: Based on the rank relationship of gene expression levels between gene pairs, gene pairs with stable expression between genes in PD samples were selected as reference stable gene pairs. Genes in PD samples were combined in pairs to form n(n-1) / 2 gene pairs. P(gene a >gene b ) ≥ 0.8 as the reference stable gene pair; (1) in m represents the total number of PD samples, k Represents the gene in PD samples a >gene b The number of samples, gene a >gene b Indicates gene a and genes b The rank relationship of expression levels; Step B2: Based on the reference stable gene pairs obtained in PD samples (gene a , geneb ), calculate the gene of each reference stable gene pair in PD samples and PRCR samples a >gene b Number of samples m 1. n 1 and gene a <gene b Number of samples m 2. n 2; Then, Fisher's exact test was used to determine whether there was a significant difference in the REO distribution of the reference stable gene pair between PD and PRCR samples, and the gene pairs that satisfied the false positive discovery rate (FDR) of less than 5% were defined as reversal gene pairs; for each reversal gene pair, the formula △P=P PD (gene a >gene b )-P PRCR (gene a >gene b ) to calculate the reversal rate △P; the larger the △P value, the greater the difference in REO of the gene pair; △P=1 means that the REO of the gene pair in PD samples is gene a >gene b , in PRCR samples, both are genes a ≤gene b ; Given a threshold △Pt such as 0.7, the gene pairs that satisfy △P>△Pt are taken as candidate prediction gene pairs.

[0008] Furthermore, the step C specifically comprises the following steps: A greedy algorithm is used to find the gene pair combination with the maximum sample coverage. Starting from each gene pair in the candidate prediction gene pair, it is used as a seed, and the remaining gene pairs are added to the current combination in turn; after each addition, the sample coverage of the new combination is calculated. This process continues until the sample coverage can no longer be increased by adding candidate gene pairs. The final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; the frequency of each characteristic gene pair in all combinations is calculated. The frequency reflects the importance of the characteristic pair, and the optimal prediction gene pair is selected as the final prediction gene pair; for a sample to be predicted, the gene in the prediction gene pair combination is calculated. a >gene b The number of gene pairs; according to the voting rule, if gene a >gene bIf the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR); then the prediction gene pair is verified using the validation set, and the obtained prediction gene pair is used as a marker for predicting the response to anti-PD-1 immunotherapy for melanoma; including: The REO reference pattern of a gene pair is defined as the gene expressed in PD samples. a >gene b Or in PRCR samples, it appears as gene a ≤gene b ; When at least one gene pair in a sample shows the gene pair REO reference pattern, the combination is considered to cover the sample, and the percentage of covered samples to the total number of samples is calculated; In order to screen a set of characteristic gene pairs that can predict the therapeutic effect, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage. Starting from each gene pair in the candidate prediction gene pair, it is used as a seed, and the remaining gene pairs are added to the current combination in turn; after each addition, the sample coverage of the new combination is calculated. This process continues until the sample coverage can no longer be increased by adding candidate gene pairs. The final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; Calculate the frequency of each characteristic gene pair in all combinations. The frequency reflects the importance of the characteristic pair. Sort the frequencies from large to small and select the ones with the highest frequency first. k feature pairs; k It is a predefined parameter with a value of 1~ n ( n is the total number of characteristic gene pairs); for the selected k feature pairs, calculate the geometric mean of their negative predictive value and positive predictive value; select the first one that reaches the maximum geometric mean and tends to be stable k gene pairs as the final prediction gene pair combination; for a sample to be predicted, calculate the gene a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR).

[0009] The present invention also relates to a method for verifying a gene pair predicting a melanoma anti-PD-1 immunotherapy response, comprising the following steps: Step A: obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from the Gene Expression Comprehensive Database and the Cancer Genome Atlas, preprocessing the expression profile data and clinical information, and dividing the melanoma expression profile data sets from different detection platforms into two data subsets, one melanoma anti-PD-1 immunotherapy data subset as a training set, and the other data subset as a validation set; Step B: Combine genes in the training set in pairs, compare the differences in relative expression rank relationships between gene pairs in PD samples and PRCR samples, screen gene pairs with significantly reversed relative expression rank relationships between genes, and identify candidate predictive gene pairs with the potential to predict the response to anti-PD-1 immunotherapy in melanoma; Step C: Use the greedy algorithm to find the gene pair combination with the maximum sample coverage; select the optimal prediction gene pair combination as the final prediction gene pair; for a sample to be predicted, calculate the gene a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR); the prediction gene pair is then verified using the validation set, and the obtained prediction gene pair is used as a marker for predicting the response to anti-PD-1 immunotherapy of melanoma; Step D: verifying the melanoma anti-PD-1 immunotherapy response prediction gene pair in the validation set, specifically including: using the melanoma anti-PD-1 immunotherapy response prediction gene pair to calculate evaluation indicators (accuracy, positive predictive value, negative predictive value, specificity, sensitivity) for the data set containing anti-PD-1 immunotherapy drug response; performing survival analysis and multivariate Cox regression analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair; performing immune infiltration analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair; and performing immunotherapy response analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair.

[0010] Furthermore, the calculation of evaluation indicators (accuracy, positive predictive value, negative predictive value, specificity, sensitivity) for a data set containing anti-PD-1 immunotherapy drug responses using melanoma anti-PD-1 immunotherapy response prediction genes described in step D specifically includes the following steps: Use a 2x2 confusion matrix to represent the prediction results of the predicted gene pairs, where the rows represent the actual labels and the columns represent the predicted labels; According to Table 1, the confusion matrix is ​​established, and the prediction index calculation formulas are: , , , ,

[0011]

[0012] Note: TP (True positive); TN (True negative); FN (False negative); FP (Falsepositive); P (positive); N (negative) Furthermore, the survival analysis and multivariate Cox regression analysis of the melanoma anti-PD-1 immunotherapy response prediction gene pairs described in step D specifically include the following steps: Kaplan-Meier (KM) survival analysis was performed using the R package "survival" and "survminer" to verify whether there was a survival difference between the two subgroups predicted by the predicted gene pairs in the validation set; the gene pairs in the predicted gene pair combination were ranked according to their expression order in the sample. a Greater than gene b The proportion of predicted gene pairs was used to construct a scoring standard for predicted gene pairs; a multivariate Cox regression analysis was performed on the validation set, incorporating factors such as age, gender, and predicted gene pair scores.

[0013] Furthermore, the immune infiltration analysis of the melanoma anti-PD-1 immunotherapy response prediction gene pair described in step D specifically includes the following steps: The ssGSEA (Single Sample Gene Set Enrichment Analysis) algorithm was used to evaluate the immune cell content between the response group and the non-response group predicted by the prediction gene pairs in the validation set. The immune cells between the two subgroups predicted by the prediction gene pairs in the validation set included 28 types: activated CD4+ / CD8+ T cells, central memory CD4+ / CD8+ T cells, effector memory CD4+ / CD8+ T cells, γδT cells, helper T cells (Th1), Th2 cells, Th17 cells, regulatory T cells (Treg), follicular helper T cells (Tfh), activated B cells, naive B cells, memory B cells, activated dendritic cells (DC), naive DC, plasmacytoid DC, natural killer cells (NK), natural killer T cells (NKT), myeloid-derived suppressor cells (MDSC), macrophages, monocytes, mast cells, eosinophils, neutrophils, CD56bright NK cells, CD56dim NK cells, Cells were compared using the Wilcoxon test.

[0014] Furthermore, the immunotherapy response analysis of the melanoma anti-PD-1 immunotherapy response prediction gene pair described in step D specifically includes the following steps: Ten published immune checkpoint-related genes (CD274 (PD-1), CTLA-4, CXCR3, LAG3, PDCD1 (PD-L1), CD80, ICOS, IFNG, IL10, and TNFRSF9) were collected, and the expression levels of these genes in the response group and non-response group predicted by the validation set were evaluated; the differences in the levels of immunotherapy response-related genes between the two subgroups divided by the prediction gene pairs in the validation set were understood by comparing the expression levels of these immune checkpoint-related genes between the two subgroups predicted by the prediction gene pairs in the validation set; and the Wilcoxon test method was used for comparison.

[0015] Compared with the prior art, the present invention has the following beneficial effects: The method for establishing and verifying a gene pair for predicting the response of melanoma to anti-PD-1 immunotherapy described in the present invention is also a method for establishing and verifying a gene pair for predicting the response of melanoma to anti-PD-1 immunotherapy. A melanoma anti-PD-1 immunotherapy data set is selected as a training set, and a gene pair matrix is ​​formed based on the rank relationship of relative expression of genes in the sample to screen candidate prediction gene pairs. Based on the candidate prediction gene pairs, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage, and the gene pair combination with the best prediction efficiency is screened as the prediction gene pair. For a sample to be predicted, the gene in the prediction gene pair combination is calculated.a >gene b According to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR). This method has been verified in the sample validation set of other different platforms. This method is based on the relative rank relationship between genes and can be robustly applied to independent clinical samples evaluated in different laboratories at an individualized level to accurately predict the occurrence and development of cancer. The predicted gene pairs can achieve truly personalized prediction of patient drug response and have clinical application value in predicting the response to anti-PD-1 immunotherapy for melanoma. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 This is a schematic diagram of a process for establishing a gene pair for predicting the response to anti-PD-1 immunotherapy of melanoma according to Example 1 of the present invention; Figure 2 This is a schematic diagram of the process of establishing a melanoma anti-PD-1 immunotherapy training set and a validation set according to Example 2 of the present invention; Figure 3 This is a schematic flow chart of a method for validating a gene pair for predicting response to anti-PD-1 immunotherapy for melanoma described in Example 2 of the present invention; Figure 4 This is a schematic diagram of the survival analysis and multivariate Cox regression analysis of genes predicting the response to anti-PD-1 immunotherapy for melanoma described in Example 2 of the present invention; Figure 5 This is a schematic diagram of immune infiltration analysis of genes predicting response to anti-PD-1 immunotherapy for melanoma according to Example 2 of the present invention; Figure 6 This is a schematic diagram of the analysis of the immunotherapy response prediction genes for melanoma anti-PD-1 immunotherapy response described in Example 2 of the present invention. DETAILED DESCRIPTION

[0017] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0018] The present application provides a method for establishing and verifying a gene pair predicting the response to anti-PD-1 immunotherapy of melanoma, which identifies biomarkers related to the response to anti-PD-1 immunotherapy based on the relative expression rank relationship between genes in a sample, and obtains a set of gene pairs as the gene pairs predicting the response to anti-PD-1 immunotherapy of melanoma.

[0019] Embodiment 1: like Figure 1 As shown, a method for establishing a gene pair for predicting the response to anti-PD-1 immunotherapy of melanoma is provided, which specifically includes: Step 102, obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from a comprehensive gene expression database and a cancer genome atlas, preprocessing the expression profile data and clinical information, and dividing the melanoma expression profile data sets from different detection platforms into two data subsets, one melanoma anti-PD-1 immunotherapy data subset as a training set, and the other data subset as a validation set; Gene expression profiles of 4 independent datasets were retrospectively collected, containing 639 samples. The analyzed data came from the Gene Expression Omnibus (GEO, http: / / www.ncbi.nlm.nih.gov / geo / ) and the Cancer Genome Atlas (TCGA, https: / / portal.gdc.cancer.gov / ). The melanoma anti-PD-1 immunotherapy expression profile dataset GSE91061 from the GPL9052 detection platform includes samples before and during treatment, which are used as training sets (GSE91061_Pre) and validation sets (GSE91061_On), respectively. The remaining three sets of melanoma expression profile datasets were used as independent validation sets, which were the melanoma expression profile data of TCGA, GSE115821, and GSE22153 of the detection platforms GPL11154, GPL18573, and GPL6102. Among them, GSE115821 also contains samples before and during treatment, while GSE22153 and TCGA-SKCM do not contain anti-PD-1 treatment information.

[0020] Before constructing the predicted gene pairs, the data needs to be preprocessed as follows: 1) Download expression profile data from the GEO database and use the platform annotation file to map the probe ID to the gene ID; 2) Remove SD samples from the training set; 3) Remove low-expression genes, which are genes whose expression is missing or zero in more than 70% of sample genes.

[0021] Step 104, combining the training set genes in pairs, comparing the relative expression rank relationship differences between the gene pairs in the PD sample and the PRCR sample, and determining candidate prediction gene pairs with potential for predicting anti-PD-1 immunotherapy for melanoma; The genes in the training set samples are combined in pairs to form gene pairs. a and gene b Representative genea and b The expression level of these two genes is a >gene b or gene a ≤gene b . Screening for genes that are maintained in more than 80% of the samples a >gene b The gene pairs were used as reference stable gene pairs for PD samples in anti-PD-1 immunotherapy for melanoma, and a total of 83,634,592 reference stable gene pairs were screened. The gene expression of each reference stable gene pair in PD samples and PRCR samples was calculated. a >gene b Number of samples m 1. n 1 and gene a <gene b Number of samples m 2. n 2. Then, Fisher's exact test was used to determine whether there was a significant difference in the REO distribution of the reference stable gene pair between PD and PRCR samples. Gene pairs with a false positive discovery rate (FDR) of less than 5% were defined as reversal gene pairs, and a total of 453,714 reversal gene pairs were screened. For each reversal gene pair, according to the formula △P=P PD (gene a >gene b )-P PRCR (gene a >gene b ) to calculate the reversal rate △P. The larger the △P value, the greater the difference in REO of the gene pair. △P=1 means that the REO of the gene pair in PD samples is gene a >gene b , in PRCR samples, both are genes a ≤gene b Given a threshold △Pt such as 0.7, the gene pairs satisfying △P>Pt are taken as candidate prediction gene pairs, and 320 reversal gene pairs are obtained as candidate prediction gene pairs.

[0022] Step 106: Based on the candidate prediction gene pairs, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage, and the gene pair combination with the best prediction performance is selected as the prediction gene pair. a >gene b According to the voting rule, if gene a >gene bIf the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR).

[0023] In order to construct a predictive gene pair for anti-PD-1 immunotherapy response in melanoma, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage based on the candidate prediction gene pairs, and the gene pair combination with the best prediction efficiency is selected as the prediction gene pair. The specific method is as follows: starting from each gene pair in the candidate prediction gene pair, use it as a seed, and add the remaining gene pairs to the current combination in turn. After each addition, the sample coverage of the new combination is calculated. This process continues until the sample coverage can no longer be increased by adding candidate gene pairs. The final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs.

[0024] Calculate the frequency of each characteristic gene pair in all combinations. The frequency reflects the importance of the characteristic pair. Sort the frequencies from large to small, and select the ones with the highest frequency first. k feature pairs. Here, k It is a predefined parameter with a value of 1~ n ( n is the total number of characteristic gene pairs). k feature pairs, calculate the geometric mean of their negative predictive value and positive predictive value. Select the first 15 gene pairs that reach the maximum geometric mean and tend to be stable as the final prediction gene pairs (15-pairs). For a sample to be predicted, calculate the gene a >gene b According to the voting rule, if gene a >gene b If the number of gene pairs of is more than half, the sample to be predicted can be classified as a responder (R), otherwise it can be classified as a non-responder (NR). The 15-pairs information is shown in Table 2.

[0025] Table 2.15-pairs

[0026] In the above-mentioned method for establishing a gene pair for predicting the response of melanoma anti-PD-1 immunotherapy, the melanoma anti-PD-1 immunotherapy data set is used as a training set, and candidate prediction gene pairs are screened based on the difference in relative expression rank relationships between gene pairs in PD samples and PRCR samples. Based on the candidate prediction gene pairs, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage, and the gene pair combination with the best prediction efficiency is selected as the prediction gene pair. For a sample to be predicted, the gene in the prediction gene pair combination is calculated. a >gene b According to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted can be classified as a responder (R), otherwise it can be classified as a non-responder (NR). This method has been verified in the validation sets of other different platforms. This method is based on the relative rank relationship between genes and can be robustly applied to independent clinical samples evaluated in different laboratories at an individualized level to accurately predict the occurrence and development of cancer. The predicted gene pairs can achieve truly personalized prediction of patient drug response and have clinical application value in predicting the response to anti-PD-1 immunotherapy for melanoma.

[0027] Embodiment 2: like Figure 2 As shown, a validation method for predicting gene pairs for anti-PD-1 immunotherapy response in melanoma is provided, specifically comprising: Step 302, obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from a comprehensive gene expression database and a cancer genome atlas, preprocessing the expression profile data and clinical information, and dividing the melanoma expression profile data sets from different detection platforms into two data subsets, one melanoma anti-PD-1 immunotherapy data subset as a training set, and the other data subset as a validation set; Step 304, combining genes in the training set in pairs, comparing the differences in relative expression rank relationships of gene pairs in PD samples and PRCR samples, screening gene pairs with significantly reversed relative expression rank relationships between genes, and determining candidate predictive gene pairs with potential ability to predict the response to anti-PD-1 immunotherapy for melanoma; Step 306: Based on the candidate prediction gene pairs, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage, and the gene pair combination with the best prediction efficiency is selected as the prediction gene pair. a >gene b According to the voting rule, if gene a >geneb If the number of gene pairs is more than half, the sample to be predicted can be classified as a responsive sample (responder, R), otherwise it can be classified as a non-responder (non-responder, NR).

[0028] Step 308, verifying the melanoma anti-PD-1 immunotherapy response prediction gene pair in a validation set, specifically including: using the melanoma anti-PD-1 immunotherapy response prediction gene pair to calculate evaluation indicators (accuracy, positive predictive value, negative predictive value, specificity, sensitivity) for a data set containing anti-PD-1 immunotherapy drug responses; performing survival analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair; performing immune infiltration analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair; and performing immunotherapy response analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair.

[0029] The evaluation indicators (accuracy, positive predictive value, negative predictive value, specificity, and sensitivity) of the data set containing the response to anti-PD-1 immunotherapy drugs were calculated using the melanoma anti-PD-1 immunotherapy response prediction gene pairs: GSE91061_Pre, GSE91061_On, GSE115821_Pre, and GSE115821_On were calculated using a 2x2 confusion matrix, and the specific verification results are shown in Table 3. The results show that 15-pairs has a good predictive effect on predicting anti-PD-1 treatment of melanoma.

[0030] Table 3. Verification results of training set and validation set

[0031] Survival analysis and multivariate Cox regression analysis were performed on the validation set using gene pairs predicting response to anti-PD-1 immunotherapy in melanoma. After removing samples with a survival time of less than 1 month, 451 samples remained in the TCGA-SKCM data set and 47 samples remained in GSE22153. In the TCGA-SKCM data set, 193 cases were predicted to be NR samples and 258 cases were predicted to be R samples. In the GSE22153 data set, 31 cases were predicted to be NR samples and 16 cases were predicted to be R samples. Figure 4 As shown in the KM survival curve, there are significant survival differences between the two types of samples in the TCGA-SKCM and GSE22153 datasets ( Figure 4A, B), and the survival of non-responding patients was worse (p<0.05, Log-rank test). This result indicates that there is a significant difference in survival time between the two subgroups predicted by 15-pairs. Multivariate Cox regression analysis was performed on GSE22153 and TCGA-SKCM, incorporating factors such as age, gender, and 15-pairs.

[0032] The results showed that 15-pairs was significantly associated with the overall survival of patients (HR=9.5, 95%CI: 3.44-26.5, P value <0.001; HR=51.1, 95%CI: 2.73-954.9, P value = 0.008, Wald test) Figure 4 C, D), the results showed that 15-pairs is an independent influencing factor for anti-PD-1 immunotherapy of melanoma. This analysis further proves the effectiveness of 15-pairs.

[0033] Based on the immune cell characteristic gene set containing characteristic genes of immune cells, the ssGSEA (Single Sample Gene Set Enrichment Analysis) algorithm was used to perform immune cell infiltration analysis on the response group and non-response group of the validation set TCGA-SKCM and GSE22153. The relative abundance of specific immune cell types in each sample was evaluated by the ssGSEA algorithm to compare the difference in the degree of immune cell infiltration between the two subgroups divided by 15-pairs. The analysis results showed that there were significant differences in T cells, B cells, dendritic cells, macrophages, eosinophils, mast cells, neutrophils, natural killer cells, etc. in the two subgroups of TCGA-SKCM divided by 15-pairs (p<0.05, Wilcoxon test), and the response group had a higher infiltration level ( Figure 5 A). At the same time, there were significant differences in B cells, dendritic cells, macrophages, eosinophils, mast cells, monocytes, myeloid-derived suppressor cells, regulatory T cells, and natural killer cells between the two subgroups of GSE22153 divided by 15-pairs (p<0.05, Wilcoxon test), and the response group had a higher infiltration level ( Figure 5 B). These results indicate that the responder groups based on 15-pairs division showed higher levels in immune cell infiltration, revealing that immunotherapy drug response may be associated with the level of immune cell infiltration.

[0034] To evaluate whether there are differences in immunotherapy between the response group and the non-response group divided by 15-pairs: ten published immune checkpoint-related genes (CD274 (PD-1), CTLA-4, CXCR3, LAG3, PDCD1 (PD-L1), CD80, ICOS, IFNG, IL10, TNFRSF9) were collected, and the expression levels of these genes in the two subgroups were evaluated. The results showed that the expression values ​​of CD274 (PD-L1), CTLA-4, CXCR3, LAG3, PDCD1 (PD-1), CD80, ICOS, IFNG, IL10, and TNFRSF9 in the response group of the TCGA-SKCM dataset were significantly higher than those in the non-response group (p<0.05, Wilcoxon test, Figure 6 A). However, since the expression of CD274 and IFNG was not detected in the GSE22153 dataset, this study analyzed the expression of the remaining 8 immune checkpoints. The results showed that the expression values ​​of CTLA-4, CXCR3, PDCD1, CD80, ICOS, IL10, and TNFRSF9 in the response group in GSE22153 were significantly higher than those in the non-response group (p<0.05, Wilcoxon test, Figure 6 B). Taking into account the prediction results of 15-pairs and the expression of immunotherapy-related genes, we can more comprehensively evaluate the differences in immunotherapy response among different subgroups, which indicates that 15-pairs can indeed predict patients' response to anti-PD-1 immunotherapy.

[0035] It should be understood that although Figure 1 and Figure 3 The steps in the flowchart are shown in sequence as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified in this document, there is no strict order restriction for the execution of these steps, and these steps can be executed in other orders. Moreover, Figure 1 and Figure 3 At least part of the steps may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least part of the sub-steps or stages of other steps.

[0036] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0037] The above embodiments only express several implementation methods of the present application, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for a person of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. Therefore, the protection scope of the patent of the present application shall be subject to the attached claims.

Claims

1. A method for establishing a gene pair for predicting response to anti-PD-1 immunotherapy for melanoma, characterized in that: The following steps are involved: Step A: obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from a comprehensive gene expression database and a cancer genome atlas, and preprocessing the expression profile data and clinical information; The melanoma expression profile datasets from different detection platforms were divided into two data subsets, one melanoma anti-PD-1 immunotherapy dataset was used as a training set, and the other data subset was used as a validation set; Step B: Combine genes in the training set in pairs, compare the differences in the relative expression rank relationships of gene pairs in disease progression PD samples and partial or complete response PRCR, screen gene pairs with significantly reversed relative expression rank relationships between genes, and identify candidate predictive gene pairs with potential ability to predict the response to anti-PD-1 immunotherapy for melanoma; Step C: Use a greedy algorithm to find the gene pair combination with the maximum sample coverage; Starting from each gene pair in the candidate prediction gene pair, use it as a seed and add the remaining gene pairs to the current combination in turn; after each addition, calculate the sample coverage of the new combination; this process continues until the sample coverage can no longer be increased by adding candidate gene pairs; the final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; Calculate the frequency of each characteristic gene pair in all combinations. The frequency reflects the importance of the characteristic pair. Select the best prediction gene pair as the final prediction gene pair. For a sample to be predicted, calculate the gene in the predicted gene pair combination a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted is classified as a responsive sample R, otherwise it is a non-responsive sample NR; the prediction gene pair is then verified using the validation set, and the obtained prediction gene pair is used as a marker for predicting the response to anti-PD-1 immunotherapy of melanoma.

2. The method for establishing a gene pair for predicting response to anti-PD-1 immunotherapy for melanoma according to claim 1, characterized in that: Acquiring expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from the gene expression comprehensive database and the cancer genome atlas as described in step A, and preprocessing the expression profile data and clinical information, specifically comprises the following steps: Step A1: Download expression profile data from the GEO database and use the platform annotation file to map the probe ID to the gene ID; Step A2: Studies have shown that in the process of anti-PD-1 immunotherapy for melanoma, the subsequent efficacy of SD samples in a stable disease state is not stable, and they will develop into PD samples or PRCR samples during subsequent treatment. Therefore, SD samples are removed from the training set; Step A3: Remove low-expression genes, which are genes whose expression is missing or zero in more than 70% of sample genes.

3. The method for establishing a gene pair predicting response to anti-PD-1 immunotherapy for melanoma according to claim 2, characterized in that: In step A1, the expression profile data is downloaded from the GEO database, and the probe ID is mapped to the gene ID using the platform annotation file. The melanoma anti-PD-1 immunotherapy expression profile dataset GSE91061 from the GPL9052 detection platform, including samples before and during treatment, is used as the training set GSE91061_Pre and the validation set GSE91061_On, respectively, and the other three sets of melanoma expression profile datasets are used as independent validation sets; the three sets of validation sets are the melanoma expression profile data of TCGA, GSE115821, and GSE22153 of the detection platforms GPL11154, GPL18573, and GPL6102.

4. The method for establishing a gene pair for predicting response to anti-PD-1 immunotherapy for melanoma according to claim 1, characterized in that: The steps described in step B are to combine the genes in the training set in pairs, compare the differences in the relative expression rank relationships of the gene pairs in the PD samples and the PRCR samples, screen the gene pairs with significantly reversed relative expression rank relationships between the genes, and determine the candidate predictive gene pairs with potential ability to predict the response to anti-PD-1 immunotherapy for melanoma, specifically including the following steps: Step B1: Based on the rank relationship of gene expression levels between gene pairs, gene pairs with stable expression between genes in PD samples were selected as reference stable gene pairs. Genes in PD samples were combined in pairs to form n(n-1) / 2 gene pairs. P(gene a >gene b ) ≥ 0.8 as the reference stable gene pair; （1） in m represents the total number of PD samples, k Represents the gene in PD samples a >gene b The number of samples, gene a >gene b Indicates gene a and genes b The rank relationship of expression levels; Step B2: Gene pairing based on reference stable genes obtained in PD samples a , gene b Calculate the gene pair of each reference stable gene in PD samples and PRCR samples a >gene b Number of samples m 1. n 1 and gene a <gene b Number of samples m 2. n 2; Then, Fisher's exact test was used to determine whether there was a significant difference in the REO distribution of the reference stable gene pair between PD and PRCR samples, and the gene pairs with a false positive discovery rate FDR less than 5% were defined as reversal gene pairs; for each reversal gene pair, the formula △P=P PD (gene a >gene b )-P PRCR (gene a >gene b ) to calculate the reversal rate △P; the larger the △P value, the greater the difference in REO of the gene pair; △P=1 means that the REO of the gene pair in PD samples is gene a >gene b , in PRCR samples, both genes a ≤gene b ; Given a threshold △Pt such as 0.7, the gene pairs that satisfy △P>△Pt are taken as candidate prediction gene pairs.

5. The method for establishing a gene pair for predicting response to anti-PD-1 immunotherapy of melanoma according to claim 1, characterized in that: The step C specifically comprises the following steps: A greedy algorithm is used to find the gene pair combination with the maximum sample coverage. Starting from each gene pair in the candidate prediction gene pair, it is used as a seed, and the remaining gene pairs are added to the current combination in turn; after each addition, the sample coverage of the new combination is calculated. This process continues until the sample coverage can no longer be increased by adding candidate gene pairs. The final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; the frequency of each characteristic gene pair in all combinations is calculated. The frequency reflects the importance of the characteristic pair, and the optimal prediction gene pair is selected as the final prediction gene pair; for a sample to be predicted, the gene in the prediction gene pair combination is calculated. a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs exceeds half, the sample to be predicted is classified as a responsive sample R, otherwise it is a non-responsive sample NR; then the prediction gene pair is verified using the validation set, and the obtained prediction gene pair is used as a marker for predicting the response to anti-PD-1 immunotherapy for melanoma; including: The REO reference pattern of a gene pair is defined as the gene expressed in PD samples. a >gene b Or in PRCR samples, it appears as gene a ≤gene b ; When at least one gene pair in a sample shows the gene pair REO reference pattern, the combination is considered to cover the sample, and the percentage of covered samples to the total number of samples is calculated; In order to screen a set of characteristic gene pairs that can predict the therapeutic effect, a greedy algorithm is used to find the gene pair combination with the maximum sample coverage. Starting from each gene pair in the candidate prediction gene pair, it is used as a seed, and the remaining gene pairs are added to the current combination in turn; after each addition, the sample coverage of the new combination is calculated. This process continues until the sample coverage can no longer be increased by adding candidate gene pairs. The final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; Calculate the frequency of each characteristic gene pair in all combinations. The frequency reflects the importance of the characteristic pair. Sort the frequencies from large to small and select the ones with the highest frequency first. k feature pairs; k It is a predefined parameter with a value of 1~ n , n is the total number of characteristic gene pairs; for the selected k feature pairs, calculate the geometric mean of their negative predictive value and positive predictive value; select the first one that reaches the maximum geometric mean and tends to be stable k gene pairs as the final prediction gene pair combination; for a sample to be predicted, calculate the gene a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted is classified as a responsive sample R, otherwise it is a non-responsive sample NR.

6. A method for validating gene pairs for predicting response to anti-PD-1 immunotherapy of melanoma, characterized in that: The following steps are involved: Step A: obtaining expression profile data and clinical information of melanoma patients and melanoma anti-PD-1 immunotherapy patients from a comprehensive gene expression database and a cancer genome atlas, and preprocessing the expression profile data and clinical information; The melanoma expression profile datasets from different detection platforms were divided into two data subsets, one melanoma anti-PD-1 immunotherapy dataset was used as a training set, and the other data subset was used as a validation set; Step B: Combine genes in the training set in pairs, compare the differences in the relative expression rank relationships of gene pairs in disease progression PD samples and partial or complete response PRCR, screen gene pairs with significantly reversed relative expression rank relationships between genes, and identify candidate predictive gene pairs with potential ability to predict the response to anti-PD-1 immunotherapy for melanoma; Step C: Use a greedy algorithm to find the gene pair combination with the maximum sample coverage; Starting from each gene pair in the candidate prediction gene pair, use it as a seed and add the remaining gene pairs to the current combination in turn; after each addition, calculate the sample coverage of the new combination; this process continues until the sample coverage can no longer be increased by adding candidate gene pairs; the final combination is the characteristic gene pair combination with the maximum sample coverage identified based on the candidate prediction gene pairs; Calculate the frequency of each characteristic gene pair in all combinations. The frequency reflects the importance of the characteristic pair. Select the best prediction gene pair as the final prediction gene pair. For a sample to be predicted, calculate the gene in the predicted gene pair combination a >gene b The number of gene pairs; according to the voting rule, if gene a >gene b If the number of gene pairs is more than half, the sample to be predicted is classified as a responsive sample R, otherwise it is a non-responsive sample NR; the prediction gene pair is then verified using the validation set, and the obtained prediction gene pair is used as a marker for predicting the response to anti-PD-1 immunotherapy for melanoma; Step D: verifying the melanoma anti-PD-1 immunotherapy response prediction gene pair in the validation set, specifically including: using the melanoma anti-PD-1 immunotherapy response prediction gene pair to evaluate the data set containing anti-PD-1 immunotherapy drug response: calculation of accuracy, positive predictive value, negative predictive value, specificity, and sensitivity; performing survival analysis and multivariate Cox regression analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair; performing immune infiltration analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair; and performing immunotherapy response analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair.

7. The method for validating a gene pair for predicting response to anti-PD-1 immunotherapy of melanoma according to claim 6, characterized in that: Step D uses melanoma anti-PD-1 immunotherapy response prediction gene pairs to evaluate the data set containing anti-PD-1 immunotherapy drug response indicators: calculation of accuracy, positive predictive value, negative predictive value, specificity, and sensitivity, specifically including the following steps: Use a 2x2 confusion matrix to represent the prediction results of the predicted gene pairs, where the rows represent the actual labels and the columns represent the predicted labels; The confusion matrix is ​​established based on the predicted genes and the prediction index calculation formulas are: 、 、 、 、 Note: TP, True positive; TN, True negative; FN, False negative; FP, False positive; P, positive; N, negative.

8. The method for validating a gene pair for predicting response to anti-PD-1 immunotherapy of melanoma according to claim 6, characterized in that: The step D of performing survival analysis and multivariate Cox regression analysis on the melanoma anti-PD-1 immunotherapy response prediction gene pair specifically comprises the following steps: Kaplan-Meier KM survival analysis was performed using the "survival" and "survminer" packages in the R package to verify whether there was a survival difference between the two subgroups predicted by the predicted gene pairs in the validation set; the expression rank relationship of the gene pairs in the predicted gene pair combination in the sample was gene a Greater than gene b The proportion of predicted gene pairs was used to construct a scoring standard for predicted gene pairs; a multivariate Cox regression analysis was performed on the validation set, incorporating factors such as age, gender, and predicted gene pair scores.

9. The method for validating a gene pair for predicting response to anti-PD-1 immunotherapy of melanoma according to claim 6, characterized in that: The immune infiltration analysis of the melanoma anti-PD-1 immunotherapy response prediction gene pair described in step D specifically includes the following steps: The ssGSEA algorithm was used to evaluate the immune cell content between the response group and the non-response group predicted by the prediction gene pairs in the validation set. The immune cells between the two subgroups predicted by the prediction gene pairs in the validation set included 28 types: activated CD4+ / CD8+ T cells, central memory CD4+ / CD8+ T cells, effector memory CD4+ / CD8+ T cells, γδT cells, helper T cells Th1, Th2 cells, Th17 cells, regulatory T cells Treg, follicular helper T cells Tfh, activated B cells, naive B cells, memory B cells, activated dendritic cells DC, naive DC, plasmacytoid DC, natural killer cells NK, natural killer T cells NKT, myeloid-derived suppressor cells MDSC, macrophages, monocytes, mast cells, eosinophils, neutrophils, CD56bright NK cells, and CD56dim NK cells; the Wilcoxon test method was used for comparison.

10. The method for validating a gene pair for predicting response to anti-PD-1 immunotherapy of melanoma according to claim 6, characterized in that: The immunotherapy response analysis of the melanoma anti-PD-1 immunotherapy response prediction gene pair described in step D specifically includes the following steps: Ten published immune checkpoint-related genes were collected: CD274 PD-1, CTLA-4, CXCR3, LAG3, PDCD1PD-L1, CD80, ICOS, IFNG, IL10, and TNFRSF9, and the expression levels of these genes in the response group and non-response group predicted by the validation set were evaluated; by comparing the expression levels of these immune checkpoint-related genes between the two subgroups predicted by the predictive gene pairs in the validation set, the differences in the levels of immunotherapy response-related genes between the two subgroups divided by the predictive gene pairs were understood; and the Wilcoxon test method was used for comparison.

Images