Development of an immune-related gene pairs index for the prognosis analysis of metastatic melanoma

Melanoma is a skin cancer with great metastatic potential, which is responsible for the major deaths in skin cancer. Although the prognosis of melanoma patients has been improved with the comprehensive treatment, for patients with metastasis, the complexity and heterogeneity of diffuse diseases make prognosis prediction and systematic treatment difficult and ineffective. Therefore, we established a novel personalized immune-related gene pairs index (IRGPI) to predict the prognosis of patients with metastatic melanoma, which was conducive to provide new insights into clinical decision-making and prognostic monitoring for metastatic melanoma. Through complex analysis and filtering, we identified 24 immune-related gene pairs to build the model and obtained the optimal cut-off value from receiver operating characteristic curves, which divided the patients into high and low immune-risk groups. Meantime, the Kaplan–Meier analysis, Cox regression analysis and subgroup analysis showed that IRGPI had excellent prognostic value. Furthermore, IRGPI was shown that was closely associated with immune system in the subsequent tumor microenvironment analysis and gene set enrichment analysis. In addition, we broken through the data processing limitations of traditional researches in different platforms through the application of gene pairs, which would provide great credibility for our model. We believe that our research would provide a new perspective for clinical decision-making and prognostic monitoring in metastatic melanoma.

and systemic treatment are difficult and not always productive in diffusion disease due to its complexity and heterogeneity 1,6,7 . Accordingly, it is urgent to find out novel personalized therapeutic strategies and biomarkers to improve and monitor the unfavourable prognosis of metastatic melanoma.
Growing studies reveals that the immune system plays a key role in the development and progression of cancers and the relationship between body immune and tumor is complicated [8][9][10] . Tumor cells express certain antigenic components, and immune cells can penetrate into tumor tissues through chemotaxis for immune clearance 8,11 . However, when tumor microenvironment is disordered, tumor cells would evade immune elimination and suppress immune response, thus leading to tumor progression 12,13 . In melanoma, existed studies have also shown that the initiation and progression of melanoma are closely related to tumor immunity 8,12 . What's more, immunotherapy in melanoma has been widely studied on this basis [14][15][16] .
ENK et al. applied interleukin-2 (IL-2) inhalation therapy to patients with pulmonary metastatic melanoma by transferring cytokines to the tumor site, which reduced the toxicity associated with systemic IL-2 administration and achieved good efficacy 17 . Interleukin-21 (IL-21) was mainly produced by T helper cell 17 (Th17) and played a key role in the development of Th17. It had strong anti-tumor activity 18 . Petrella et al. proved the activity of IL-21 in metastatic melanoma through a multicenter phase II study, and it had a certain curative effect on metastatic melanoma 19 . High expression of cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) was closely related to antigen-specific T cell dysfunction in metastatic melanoma. Relevant preclinical studies showed that the introduction of inhibitory antibodies to CTLA-4 could eliminate downstream inhibitory signals, thus avoiding the dysfunction of antigen-specific T cell, which could produce a cytotoxic anti-tumor response 20,21 . Meanwhile, a clinical trial in the antitumor activity of programmed cell death ligand 1/programmed cell death 1(PD-L1/ PD-1) signaling blocking was confirmed useful in multiple types of cancers, including advanced melanoma 22 . These studies proved that tumor immunity was closely related to the progression and treatment of metastatic melanoma, so a new personalized comprehensive prognosis index based on immunity would be very promising.
In our study, we attempt to establish a novel prognosis index to predict the prognosis of patients with metastatic melanoma, which is based on the screening of gene pairs and could greatly reduce the biological heterogeneity and technical bias of different cross-sequencing platforms. We hope to provide a novel insight for prognosis prediction and clinical decision-making for metastatic melanoma patients.

Materials and methods
Data acquisition. The flowchart of IRGPI establishment and validation is presented in Fig. 1. The fragments per kilobase of exon model per million mapped reads (FPKM) RNA-seq data and clinical data of melanoma samples were obtained from The Cancer Genome Atlas (TCGA) data portal (https ://cance rgeno me.nih. gov/). The samples with no metastatic foci, no follow-up information, or follow-up time less than 30 days were excluded. Then, 354 metastatic melanoma samples were reserved for training dataset. Additionally, the expression profiles of GSE65904 were downloaded for testing dataset through GEOquery R package. The same sample filtering standard was used for GSE65904 dataset and 186 metastatic melanoma samples were retained. Finally, 1811 unique immune-related genes (IRGs) were acquired from ImmPort database (https ://immpo rt.niaid .nih. gov) to construct immune-related prognostic signature. Data preprocessing. The Ensembl IDs of RNA-FPKM data were transformed into gene symbols based on the Ensembl database (http://asia.ensem bl.org/index .html). The probe IDs of GSE65904 expression profiles were converted into gene symbols through the illuminaHumanv4.db R package. Ensembl IDs or probe IDs were retained on basis of the mean overall expression of each gene. Next, only IRGs measured by all platforms with relatively high variation (determined by MAD > 0.5, MAD: median absolute deviation) were selected for further analysis.
Establishment of immune-related gene pairs index for prognosis prediction. The immunerelated gene pairs (IRGPs) were constructed by pairwise comparison of the gene expression level in a specific sample or profile. If the expression level of the first IRG was higher than that of the second IRG, the score of this IRGP was 1; otherwise, the score was 0. Next, the establishment of immune-related gene pairs index (IRGPI) was based on previous description 23 . The RNA-FPKM data was identified as training dataset for establishing the IRGPI through using Lasso Cox proportion hazards regression with tenfold cross validation (glmnet R package, version: 3.0-2). The IRGPs of training dataset with a small variation and imbalanced distribution (MAD = 0) were excluded from the analysis. Then, we conducted 1000 times Lasso Cox proportion hazards regression analysis into training dataset, in which the model with the most occurrences was identified as the most stable gene pairs model and was used for the development of IRGPI.
Kaplan-Meier curve analysis and validation of IRGPI. The prognosis risk of metastatic melanoma patients was distinguished based on IRGPI and the time-dependent receiver operating characteristic (ROC) curve analysis within five years was performed to determine the optimal cut off in the training dataset. KM curve analysis and subgroup analysis were further used to evaluate the ability of IRGPI to distinguish survival risk of metastatic melanoma. Univariate and multivariate Cox analyses were conducted to compare the survival impact of IRGPI with other clinical characteristics. Further, the prognosis ability of IRGPI was validated in the independent GSE65904 testing cohort by KM curve analysis.
Tumor microenvironment in different IRGPI risk groups. The  www.nature.com/scientificreports/ infiltrating score in each sample by using 29 immune gene sets 24 . The immune score of each sample was calculated by ESTIMATE R package 9 . We used the t test to identify the difference of tumor immune microenvironment between the high and low-IRGPI risk groups.
Immune related biological processes in different IRGPI risk groups. The ordered gene lists of two cohorts were identified through corresponding R packages (training cohort: edgeR R package, testing cohort: limma R package) and Gene Set Enrichment Analysis (GSEA) 25 was conducted on the gene with false discovery rate (FDR) less than 0.05 to identify biological processes that were differently activated between the high and low-IRGPI risk groups. The number of random sample permutations were set at 1000 and the min size of gene set was set at 160. The biological processes with FDR < 0.01 in both cohorts indicated a statistically significant difference.
Statistical analysis. All statistical analyses were conducted using the R software (version 3.6.3) (http:// www.r-proje ct.org/) and its corresponding R packages. The Area Under Curve (AUC) of ROC curve was performed using the survialROC R package. KM curve analysis was completed using log-rank test from survminer R package. GSEA analysis was completed using clusterProfiler R package. The c5.bp.v7.1.entrez.gmt file from the Molecular Signatures Database (MSigDB, http://softw are.broad insti tute.org /gsea/index.jsp) was obtained for the GSEA to identify biological processes. All P value of less than 0.05 indicated a statistically significant difference in all analysis.

Results
Establishment of immune-related gene pair index for prognosis prediction. In the analysis, 376 IRGs of 354 TCGA metastatic melanoma samples were retained for constructing the IRGPs. The detail clinical features of TCGA cohorts were shown in Table 1. After removing IRGPs with relatively small variation (MAD = 0), 119 IRGPs were left for initial candidate IRGPs. Next, 24 IRGPs with the highest frequency were www.nature.com/scientificreports/ selected for the development of IRGPI through the 1000 times Lasso Cox proportion hazards regression in training cohort ( Fig. 2A). The information of 24 IRGPs and 45 unique IRGs from IRGPs was shown in Table 2. . The subgroup analysis also showed that the prediction impact of IRGPI on OS of metastatic melanoma was promising in training cohort ( Table 3). The IRGPI was also performed into the independent GSE65904 testing cohort to validate its accuracy and prediction ability. The patients of testing cohort were divided into high and low-IRGPI risk group based on the same cut off. The KM curve analysis indicated that patients with the high-IRGPI risk also had poor OS compared with the low-IRGPI risk group (Fig. 4 Different IRGPI risk groups displayed differential tumor microenvironment. The relative scores of 29 immune cells for each patient were estimated by ssGSEA algorithm. The comparative summary of ssGSEA output result was performed in these two risk groups and a wide variety of differential immune infiltration cells existed. Moreover, both the training cohort and the testing cohort showed the low-IRGPI risk group had a better total immune score than the high-IRGPI risk group by ESTIMATE algorithm (Fig. 5B , p = 2.22e−16; Fig. 5D , p = 6.2e−10).

Immune related biological processes in different IRGPI risk groups. In GSEA analysis, various
immune-associated biological processes were enriched, such as adaptive immune response; immune system development; leukocyte mediated immunity; positive regulation of immune response; activation of immune response; cell activation involved in immune response; immune response regulating signaling pathway; immune response regulating cell surface receptor signaling pathway (Fig. 6).

Discussion
With the development of researches on tumor immunity microenvironment, immunotherapy has been listed as a successful treatment option for a wide variety of cancers, including metastatic melanoma 26,27 . As an emerging and effective therapeutic approach, tumor immunotherapy has a broad prospect in metastatic melanoma 14,16,28 . In our study, we defined 376 IRGs and 119 IRGPs that were closely related to metastatic melanoma, which would provide us with powerful conditions to establish a novel prognosis model for patients with metastatic melanoma based on immunogenomic landscape analysis.
Through further analysis, 24 IRGPs, including 45 unique IRGs, were obtained, and the IRGPI of metastatic melanoma was further established based on the 24 IRGPs. The optimal cut off for distinguishing patients into high and low-IRGPI risk groups was identified based on the time-dependent ROC curve within five years. Further, KM curve analysis was used to examined the effect of IRGPI on the prognosis prediction of metastatic melanoma, and the results showed that the patients with high-IRGPI risk were associated with worse OS and DFS, www.nature.com/scientificreports/    www.nature.com/scientificreports/ indicating a poor prognosis in high-IRGPI risk group. Meanwhile, we also proved that IRGPI could be used as a promising and independent prognostic model by multivariate Cox analysis and subgroup analysis. Moreover, we used immune gene pairs for data analysis, which would not consider the technical deviation of different platforms to the greatest extent, and successfully solved the problem of different data platforms for expression. This would bring great hope for more accurate prognosis prediction in metastatic melanoma. Therefore, it is evident that IRGPI has a promising potential to be employed as a reliable prognostic index for metastatic melanoma. The IRGPI prognosis index consisted of 45 unique IRGs, and many of them have been shown to be strongly associated with cancers development. The latest research reported that IRF-1 could be used as an indicator of PD-L1 expression capability in anti-PD-1 therapy, so it could predict the therapeutic effect of metastatic melanoma 29 . Similarly, IDO-1 was also reported that was closely related to the anti-PD-1 response of metastatic Figure 5. Tumor immune microenvironment status within differential immune risk groups by using ssGSEA. (A) Summary of the 29 immune cells score of differential immune risk groups in training cohort. (B) The total immune score of differential immune risk groups in training cohort. (C) Summary of the 29 immune cells score of differential immune risk groups in testing cohort. (D) The total immune score of differential immune risk groups in testing cohort. All p values were estimated by t test (*p < 0.05, **p < 0.01, ***p < 0.001). Wente et al. proved that the expression of CXCL14 in pancreatic cancer was significantly higher than normal pancreatic tissue, especially in the frontier tissue of invasive pancreatic cancer, which indicated that CXCL14 could play an important role in tumor metastasis 32 . Additionally, related researches showed the over-expression of ITGAV was closely related to the development of colorectal cancer and spreading of colorectal cancer cells via perineural invasion 33 . And the high expression of PIK3CD was also proved that affected the distant metastasis and poor prognosis of colon cancer 34 . These studies indicated that the IRGs in IRGPI played an important role in the occurrence and development of cancers, suggesting that our prognostic model would have good prospect and prognostic value in metastatic melanoma. The infiltration of immune cells in the tumor microenvironment played an important role in the progression of cancers, and it was also reported in metastatic melanoma in a wide variety studies 27,35 . Therefore, we calculated the score of 29 immune cells infiltration in each sample using ssGSEA to further explore the correlation between tumor microenvironment immune cells infiltration and IRGPI. The results showed that different IRGPI risk groups displayed differentiates tumor immunity microenvironment. The immune infiltration cells www.nature.com/scientificreports/ played different roles in metastatic melanoma. For example, CD8+TIL has been reported for the treatment of refractory metastatic melanoma patients and obtained approved curative effect 36 . Romero et al. found that the increase of CD4+NKG2D+Th1 in patients with metastatic melanoma was associated with prolonged survival 37 .

Scientific Reports
Another study suggested that IL-2 treatment could restore Th1 ⁄ Th2 balance in metastatic melanoma and activate lymphocytes. At the same time, it enhanced the ability of monocytes producing IFN-γ, which induced a systemic immune response, thus obtained improved prognosis and better clinical benefits 38 . This is consistent with our results. Furthermore, in the pulmonary metastasis of melanoma in mice, Saga et al. proved that local generated melanoma-specific CTLs could significantly reduce the number of metastatic foci 39 . These studies confirmed that immune cells differential infiltration could be a potential mechanism of IRGPI prognosis prediction role, but these observations could be further explored to fully understand the subtle differences in microenvironment immune cell infiltration. The results strengthened our understanding of immune cells infiltration in the metastatic melanoma tumor microenvironment, which provided the conditions for further exploration.
We conducted GSEA to explore the biological processes of IRGPI acting on metastatic melanoma and demonstrated that large number of immune-related processes were differentially enriched between high and low IRGPI risk groups. Existing studies have shown that positive regulation of the immune response could reduce the risk of death from lymphatic metastasis of melanoma, and it was closely related to the increase and activation of tumor infiltrating lymphocytes (TIL) 40 . Immune response has shown to be involved in tumor development by modulating cell surface receptor signaling pathways. For example, the dynamic variations of the cell surface receptor C-X-C Motif Chemokine Receptor 3 (CXCR3) played an important role in the metastasis of melanoma 41 . Leukocyte mediated immunity could inhibit the metastasis of melanoma, such as cytotoxic T cell (CTL), which could reduce the metastatic foci 39 . This was consistent with the previous understanding. In addition, the immune response regulating signaling pathway may play a role in tumor progression and metastasis through the regulation of PD-L1 expression 42 . These results suggested that immune-related signaling pathways were an important aspect in affecting tumor progression. Our study elucidated the relevance of IRGPI in multiple immune biological processes, providing theoretical support for the application of IRGPI and the immune-related studies of subsequent metastatic melanoma.

Conclusion
The study breaks the limitations of traditional studies and reduces the biological heterogeneity and technical bias of different cross-sequencing platforms through the application of the gene pairs, which provided great convenience and reliability for the establishment of the prognostic index. Although our research showed advantages in data processing, there were still some limitations to consider. First, due to the lack of in vitro or in vivo experiments, the reliability of our molecular mechanism analysis results was limited. Second, this study was a retrospective study and a well-designed clinical trial was needed to further verify our results, though the mechanism of action of IRGPI in metastatic melanoma has been effectively elucidated through multiple methods. Therefore, the findings would provide new insights into the clinical decision-making and prognostic monitoring of metastatic melanoma and provide theoretical support for further researches. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.