From March 2016 to March 2019, 55 melanoma patients treated with anti-PD-1 therapy were recruited for this study from PUCH. Formalin-fixed, paraffin-embedded pretreatment tumor samples were obtained from all patients. We separated all the clinical and pathological data by medical record review, including sex, age, primary site, metastasis status, and clinical efficacy. Tumor responses were evaluated using the Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1, including complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). CR and PR were regarded as responders, while PD and SD were regarded as non-responders. In this study, overall survival (OS) and progression-free survival (PFS) were used as the primary and secondary survival endpoints, respectively. Gene expression data for the PUCH cohort was based on the Illumina NovaSeq 6000 platform. The details of processing the GEP of the PUCH cohort have been described in our previous study (30). This study was conducted according to the Declaration of Helsinki Principles and approved by the Medical Ethics Committee of PUCH. Informed consent for the use of material in medical research was obtained from all participants.

We obtained RNA-seq and clinical data from four publicly available cohorts of melanoma patients treated with anti-PD-1 or anti-CTLA-4 therapy, including Gide19 (n = 41) (31), Hugo16 (n = 28) (32), Riaz17 (n = 51) (33), and VanAllen15 (n = 42) (34). Data of Gide19 cohort (PRJEB23709) were downloaded from the European Nucleotide Archive database (ENA; https://www.ebi.ac.uk/ena). Data of Hugo16 cohort (GSE91061) and Riaz17 cohort (GSE78220) were downloaded from the Gene Expression Omnibus database (GEO; https://www.ncbi.nlm.nih.gov/geo/). Data of VanAllen15 cohort (phs000452.v2.p1) were downloaded from the database of Genotypes and Phenotypes (dbGap; http://www.ncbi.nlm.nih.gov/gap). The treatment response to immunotherapy consisted of CR/PR/SD/PD according to RECIST 1.1 guidelines, which were used in our analysis.

Moreover, we downloaded the RNA-sequencing data of all available cutaneous melanoma samples from The Cancer Genome Atlas (TCGA) database through the GDC tool (http://portal.gdc.cancer.gov/). The survival data of these patients were extracted from cbioPortal (http://www.cbioportal.org). Patients with OS less than one month were excluded from our analysis. In addition, we separated the TMB data of melanoma patients in the TCGA-SKCM cohort from The Cancer Immunome Atlas (http://tcia.at/home) (35).

We constructed a predictive signature based on IRGs gathered from the ImmPort Web portal (https://www.immport.org/home) (36). Two IRGs constituted one IRGP and formed as “IRG-A|IRG-B”. The score of IRGPI was generated through pairwise comparison of gene expression levels in specific samples. When the expression level of IRG-A was higher than IRG-B, the IRGP was assigned a score of 1; otherwise, the IRGP score was 0. IRGPs with score of 0 or 1 in over 80% of the specimens were regarded as IRGPs with constant values, which does not contribute to the difference of patient survival (37). Therefore, we excluded these IRGPs with constant values from our analysis.

The Gide19, Hugo16, and Riaz17 cohorts merged into the meta cohort, which was used for the construction of IRGPI. Firstly, we used the log-rank test to investigate the correlation of each IRGP to patients’ OS in the meta cohort. According to the analysis results, IRGPs with a false discovery rate (FDR) < 0.001 were candidates to build the IRGPI. Then, the multivariate Cox regression analyses were performed to obtain the hub IRGPs and the respective coefficients. Finally, the IRGPI formula was defined as follows:

The predictive and prognostic values of IRGPI for immunotherapy were validated in PUCH and VanAllen15 cohorts. Based on the treatment response to immunotherapy, we conducted the receiver operating characteristic (ROC) curve analysis to estimate the prediction accuracy of IRGPI. Using the cut-off value that generated the maximum Youden index (38), the patients were divided into IRGPI-high and IRGPI-low groups. Then, the log-rank tests were conducted for comparison of the survival outcomes between two IRGPI groups.

The transcriptomic data of TCGA-SKCM (skin cutaneous melanoma) cohort were used for analyzing the association of IRGPI with immune-related features. Using the IRGPI formula, we calculated the IRGPI score of each patient in the TCGA-SKCM cohort. The cut-off value for the IRGPI was determined on the basis of the association with patients’ OS by using X-tile software (version 3.6.1) (39). Based on two bioinformatic analyses of GEP data in the TCGA-SKCM cohort, we calculated the enrichment of immune cells between two IRGPI groups. Briefly, we used Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) method to calculate the immune score and ESTIMATE score of patients (40). CIBERSORT was further used to distinguish 22 immune cells, such as T cell types, B cell types, NK cells, and myeloid cell types (41). In the signaling analysis, we conducted gene set enrichment analysis (GSEA) to distinguish which immune-related pathways were markedly different between IRGPI-high and IRGPI-low groups.

To further characterize the tumor immune microenvironment between two IRGPI groups, we performed single simple GSEA on some previously published immune-related signatures (19, 42–45) and compared the score between IRGPI-high and IRGPI-low groups.

All statistical analyses were performed using the R software (version 3.6.3) and Prism 8. Survival analyses were performed using the R packages “survival” and “survminer”. The signature of IRGPs was obtained using the R package “glmnet”. Univariate and multivariate Cox regression analyses were conducted using the R package “survival”. ROC curve analyses were performed using the R package “survivalROC”. ESTIMATE analysis was conducted using the R package “estimate”. CIBERSORT analysis was processed using the R packages “e1701”, “preprocessCore”, and “limma”. All statistical analyses were two-sided, and P < 0.05 was considered as statistically significant.

The flowchart of this study design is presented in Figure 1 . A total of 217 patients treated with ICI from five cohorts were included in this study. We constructed the IRGPI based on the meta cohort (n = 120), which consisted Gide19 (n = 41), Hugo16 (n = 28), and Riaz17 (n = 51) cohorts. The PUCH (n = 55) and VanAllen15 (n = 42) cohorts were used for validation the prediction model. The clinicopathological characteristics are summarized in Table 1 . The median follow-up is 13.1~32.7 months in five cohorts. Notably, 43.6% of patients in PUCH cohort were acral melanomas, which have been reported to be the main subtype of melanoma in Asians. However, the vast majority of patients in other cohorts were cutaneous melanoma.

Among the 2487 IRGs from the ImmPort database, 1138 IRGs commonly occurred in the GEP of five cohorts and 1295044 IRGPs were calculated. We excluded 1251102 IRGPs (96.6%) with constant values in any cohort and 43942 IRGPs remained for subsequent construction of the IRGPI ( Table S1 ). According to the univariate Cox regression analyses of the correlation between each IRGP and patients’ OS in the meta cohort, 236 IRGPs with adjusted P < 0.001 were selected as prognostic IRGPs ( Table S2 ). Then, we performed the multivariate Cox regression analyses to filtrate IRGPs to construct the IRGPI. Finally, 31 IRGPs were filtrated to define the IRGPI ( Table 2 and Figure S1A ). The IRGPI consisted of 51 unique IRGs, most of which encoded molecules involved in antimicrobials, cytokines, and cytokine receptors.

Based on the IRGPI formula, we calculated each patient’s IRGPI score in the meta cohort and exhibited the result in a heatmap ( Figure 2A ). We conducted ROC curve analysis to evaluate the prediction accuracy of IRGPI for the efficacy of ICI therapy in the meta cohort. With a cut-off value based on the Youden index of -1.221, we found that IRGPI-low group patients were correctly classified as CR/PR with a sensitivity of 71.7% (38/53). Further, IRGPI-high group patients were successfully classified as SD/PD with a specificity of 91.0% (61/67). The overall accuracy of IRGPI was 82.5% (99/120) with AUC of 0.854 ( Figure 2B ). IRGPI-low group showed a higher ORR than IRGPI-high group (71.7% vs. 9.0%; Figures 2C, E ). Moreover, we evaluated the relationship between the IRGPI score and OS/PFS in the meta cohort. Kaplan-Meier survival analysis revealed that IRGPI-low group patients had significantly longer OS (P < 0.001; Figure 2D ). The median PFS for IRGPI-high group patients was markedly shorter than IRGPI-low group patients in the Gide19 cohort (P < 0.001; Figure S2A ). The pooled hazard ratio (HR) along with 95% confidence interval (CI) for the association between high IRGPI score and OS in 119 cases of patients was 8.02 (3.91~15.19), and no significant heterogeneity among the three datasets was observed (I² = 0%, P = 0.98, Figure S1B ). Overall, the IRGPI showed a superior prediction for the benefit of ICI therapy in the meta cohort.

To verify the robustness of IRGPI in predicting the efficacy of ICI therapy, we assessed the correlation of IRGPI score with overall response rate and survival outcomes in VanAllen15 and PUCH cohorts. In the VanAllen15 cohort, the IRGPI successfully identified 31 of 42 patients with an overall accuracy of 73.8% and an AUC of 0.767 ( Figure 3A ). Similarly, in the PUCH cohort, the IRGPI demonstrated an overall accuracy of 72.7% (40/55) and AUC of 0.737 ( Figure 3B ). The ORR of IRGPI-high group was lower than IRGPI-low group in VanAllen15 cohort (64.3% vs. 7.1%; Figures 3C and S3A ) and PUCH cohort (47.4% vs. 13.9%; Figures 3D and S3B ), respectively. As expected, higher IRGPI conferred poor survival outcomes in VanAllen15 cohort (OS: P < 0.001, PFS: P < 0.001; Figures 3E and S2B ) and PUCH cohort (OS: P = 0.004, PFS: P = 0.015; Figures 3F and S2C ), respectively. These results confirmed that the IRGPI is reliable for the prediction of ICI therapy responsiveness in VanAllen15 and PUCH cohorts.

Reportedly, the infiltration of immune cells, especially CD8⁺ T cells, is associated with immunotherapy response in many types of cancer (46). Based on the above results, we further investigated the relationship between IRGPI and tumor immune microenvironment features in melanoma patients. The TCGA-SKCM cohort was stratified into IRGPI-high and IRGPI-low groups using X-tile software ( Table S3 and Figure S4 ). Firstly, we calculated the immune score and ESTIMATE score of patients in TCGA-SKCM cohort by ESTIMATE algorithm. The data showed that both the immune score and ESTIMATE score were considerably increased in IRGPI-low group compared with IRGPI-high group ( Figures 4A, B ). Secondly, the CIBERSORT analytical tool was used to estimate the proportions of 22 types of immune cells in each SKCM sample. The results revealed that the infiltration levels of CD8⁺ T cells, activated memory CD4⁺ T cells, naive B cells, and NK cells in IRGPI-high group were lower than that in IRGPI-low group, while resting memory CD4⁺ T cells, M0 and M2 macrophages showed the opposite trend ( Figure 4C ). Finally, we performed GSEA to identify which pathways were enriched at specific IRGPI levels. As shown in Figure 4D , the pathways of inflammatory response, interferon response, antigen processing and presentation, and T cell receptor signaling were markedly upregulated in IRGPI-low group.

A series of potential biomarkers have been developed to predict the response of ICI therapy in malignant tumors, such as TMB, immune inhibitory receptor expression levels. We analyzed the relationship between IRGPI and TMB in TCGA-SKCM cohort and the results showed that higher IRGPI conferred lower TMB ( Figure 5A ). As expected, the expression levels of immune inhibitory receptors (including PD-1, CTLA-4, LAG3, TIM-3, and TIGIT) showed the same trend as TMB between IRGPI-high and IRGPI-low groups ( Figure 5B ). Moreover, the deficiency of human leukocyte antigen (HLA) could impair antigen presentation and initiate antitumor immunity, which consequently resulting in primary resistance to immunotherapy (47). We then investigated the correlation of IRGPI to the expression levels of HLA members and the data indicated most HLA members were substantially upregulated in IRGPI-low group compared with IRGPI-high group ( Figure 5C ).

Some immune-related GEP signatures have been described to predict the benefit of ICI therapy in melanoma ( Table S4 ). We therefore compared these GEP signature scores between IRGPI-high and IRGPI-low groups. Consistent with other biomarkers, these GEP signature scores were significantly downregulated in IRGPI-high group compared with IRGPI-low group ( Figure 5D ).