Identification of immune-related genes and integrated analysis of immune-cell infiltration in melanoma

Objective: This study was conducted to screen out immune-related genes in connection with the prognosis of melanoma, construct a prognosis model and explore the relevant mechanisms. Methods and materials: 1973 genes associated with immune system were derived from the Immport database, and RNA-seq data of melanoma and information of patients were searched from the Xena database. Cox univariate analysis, Lasso analysis and Cox multivariate analysis were used to screen out six genes to construct the model. Then the risk scores were estimated for patients based on our constructed prognosis model. Estimate was used to affirm that the model was about immune infiltration, and CIBERSORT was used to screen out immune cells associated with prognosis. TIDE was applied to predict the efficacy of immunotherapy. Finally, GSE65904 and GSE19234 were used to confirm the effectiveness of the model. Results: ADCYAP1R1, GPI, NTS might cause poor prognosis while IFITM1, KIR2DL4, LIF were more likely conductive to prognosis of melanoma patients and a model of prognosis was constructed on the basis of these six genes. The effectiveness of the model has been proven by the ROC curve, and the miRNAs targeting the screened genes were found out, suggesting that the immune system might impact on the prognosis of melanoma by T cell CD8+, T cell CD4+ memory and NK cells. Conclusions: In this study, the screened six genes were associated with the prognosis of melanoma, which was conductive to clinical prognostic prediction and individualized treatment strategy.


INTRODUCTION
Melanoma, a prevalent cutaneous malignancy, arises from melanocytes and exhibits a propensity for metastasis to diverse organs via lymphatic pathways [1,2].Heightened aggressiveness has been identified as a primary contributor to mortality in skin cancer cases [3,4].Middle-aged individuals are predisposed to melanoma, with males exhibiting a higher incidence than females [5,6].Recent years have witnessed an overall rise in melanoma incidence, with varying growth rates observed across different age cohorts [6].Previously, the primary risk factor for the development of tumors was believed to be excessive exposure to ultraviolet rays [7].Subsequent research has revealed that the presence of melanocytic nevi, family history, and genetic susceptibility are also significant risk factors [8].Tumorigenesis is a multifactorial process involving the intricate interplay of various factors, including genes, epigenetics, and the environment [9].The primary treatment approach for melanoma entails wide local excision, which may be complemented by chemotherapy and targeted therapy.Only a few areas that could not be operated on are treated with radiation therapy due to the radio resistance of melanoma [10].
The immune system assumes a pivotal role in the initiation of tumor formation [11].During the early stages, melanoma impedes the elimination of cancer cells through two mechanisms of tumor immune evasion.AGING Subsequently, in the advanced stages, cancer cells detach from the primary site via tumor invasion and metastasize to distant regions of the body [12,13], thereby exacerbating the disease and influencing the prognosis.Tumor infiltrating lymphocytes serve as a crucial role in the immune system's defense against tumor metastasis [14].Melanoma exhibits a high frequency of genetic mutations [15,16], resulting in the generation of neoantigens that trigger the host's immune response against the tumor.These neoantigens can be utilized for diagnostic, therapeutic, and prognostic purposes [17].Consequently, personalized care and treatment strategies can be tailored to enhance patients' physical and psychological well-being and facilitate their recovery.
This study conducted an analysis of immune genes linked to melanoma, resulting in the identification of six genes that were utilized to develop a prognostic model aimed at personalized treatment and prognosis enhancement.The accuracy of this prognostic model was also verified.In addition, the analysis has also revealed that the immune system might regulate the prognosis of melanoma through 'the regulation of lymphocyte activation', 'activation of T cell', 'cytotoxicity of NK cell', and 'differentiation of Th1 and Th2 cell' pathways.

Data and patients
1793 genes associated with immunity were extracted from the Immunology Database and Analysis Portal (Immport, http://www.immport.org/home).The RNAsequencing (RNA-seq) data of melanoma Fragments Per Kilobase Million (FPKM) in the Cancer Genome Atlas (TCGA) was obtained from the Xena database [18] and the amount of 455 melanoma patients with complete survival data was screened based on the relative information of patients.

Data analysis
According to patients' gene expression level and survival information, we identified genes that were relevant to prognosis of melanoma from 1793 immune-related genes via Cox univariate analysis, and then we establish a model which could presage the patients' prognosis of melanoma by combining the Least Absolute Shrinkage and Selection Operator (LASSO) analysis [19,20] and Cox multivariate analysis (Supplementary Material 1).

Differential expression analysis
Depending on the prognostic model, risk scores were performed for the patients, and we divided 455 people with melanoma into high risk score (RS) and low RS groups stood on the median of risk scores.The analysis of differential expression was executed on the two groups through the "LIMMA" package [21], and enrichment analysis of differentially expressed genes was executed by the "Clusterprofiler" package [22,23] (|logFC≥1, P<0.05|).
The miRNA data set of TCGA's melanoma was downloaded from Xena database, and people with melanoma were split into two groups following the above groups.Differential expression analysis was performed to construct the ceRNA network according to the principle of the ceRNA hypothesis [24] and the miRWalk [25,26] database (|logFC≥0.7,P<0.05|).
The gene mutation data set of TCGA's melanoma was downloaded from Xena database, and samples were split into two groups by risk scores.The different of mutation between high RS group and low RS group was analyzed by the "maftools" package [27].

Immune infiltration and treatment
Estimation of Stromal and Immune cells in malignant tumor tissues using Expression data (ESTIMATE) is a biological information software.In the cancer microenvironment, immune cells and stromal cells are the main normal cells in tumor tissues, and ESTIMATE can use gene expression characteristics to speculate the contents of tumor cells and different infiltrations of normal cells [28].
CIBERSORT (https://cibersort.stanford.edu/) is a method on the ground of linear support vector regression, which can describe the abundance of different cell subsets in complicated tissues from gene expression profiles [31].
Tumor Immune Dysfunction and Exclusion (TIDE) is the tool to estimate the potential of tumor immunologic escape from the gene expression profiles of cancer samples.The score of TIDE computed for each sample can predict response to immune checkpoint blockade as a biomarker.
Patients in the high RS group and the low RS group were evaluated by the ESTIMATE score, and their stromal cell immune cells and tumor purity were AGING obtained in the light of their gene expression.And we explored the relationship between risk score and TME.Then the CIBERSORT analysis was fulfilled on the two groups to compare 22 kinds of immune cells' infiltration.The analysis of TIDE was used to predict the efficacy of immunotherapy in the two groups.Statistical significance was examined using Students t -tests.

Validation
The GEO database GSE65904 and GSE19234 were chosen as validation data sets [13,32].Then people with melanoma were divided into two groups depended on the model of prognosis, and survival analysis and ROC test on the high RS group and the low RS group were performed to verify the effectiveness of the model.

Establishment of prognostic model
Cox univariate regression analysis was availed to select genes corresponding to the prognosis of melanoma from 1793 immune-related genes, while Lasso regression analysis to further curtail the number of genes.And then, Cox multivariate regression analysis was used for analyzing the outcome of melanoma.Six genes closely associated with the prognosis of melanoma patients, including ADCYAP1R1, GPI, IFITM1, KIR2DL4, LIF and NTS, were elected to fabricate a prognosis model.The calculation equation of the risk score is as shown below: risk score= (-0.30899*KIR2DL4) + (-0.13892*IFITM1) + (0.28977*GPI) + (-0.16050*LIF) + (0.48867* ADCYAP1R1) + (0.11124*NTS).The consistency index (CI) of this model was 0.67 (Figure 1A).The survival curves of these six biomarkers were consistent with the positive and negative coefficients in the formula (Figure 1B-1G).

Risk score and ROC curve analysis
Each patient received their risk score based on the above model, while their survival situation was marked, and the expression of each gene in the model from patients was shown by a heat map (Figure 2A-2C).Obviously, ADCYAP1R1, GPI and NTS might lead to poor prognosis for melanoma patients, while IFITM1, KIR2DL4 and LIF are more likely to improve outcomes.At the same time, we also analyzed the effectiveness of the model by the ROC time-dependent curve, within the Area Under the Curve (AUC) of 12, 24 and 60 months scored of 0.69, 0.72 and 0.69, respectively (Figure 2D).

Enrichment analysis
The risk score which each patient with melanoma received respectively was sorted from smallest to largest, and in consonance with the median, the patients were separated into two groups, a high RS group and a low RS group.14 genes highly expressed in the high RS group and 249 genes in the low RS group were chalked up through analysis of gene differential expression (| log FC |≥1, P<0.05).Enrichment analysis of these genes have shown that the pathways are mainly enriched in 'the regulation of lymphocyte activation', 'activation of T cell', 'cytotoxicity of NK cell', and 'differentiation of Th1 and Th2 cell' pathways (Figure 3).

CeRNA network and somatic mutations
Hinged on the above risk score grouping, the samples from the miRNA dataset were diverged into the high RS and low RS group for differential expression analysis.20 miRNAs strongly expressed in the high RS group and 18 miRNAs in the low RS group were screened by analysis of gene differential expression (| log FC |≥0.7,P<0.05) (Figure 4A, 4B).According to the theory of ceRNA hypothesis and miRWalk, targeted relationships between miRNAs with differential expression and these 6 biomarkers were found (Figure 4C).
In order to identify the difference of somatic mutations between the two sets, mutation data were envisaged by the "maftools" package in R software [27].The Variant Allele Frequencies of top mutated genes from the high RS group were stronger than those in the low RS group, suggesting the low-risk patients have better prognosis (Figure 5A, 5B).And the forest plot showed specific mutant genes in two groups (Figure 5C).

Immune infiltration analysis
In the above two groups which were separated by ESTIMATE analysis, it was found that the low RS group might have a better Immune score, Stromal score and ESTIMATE score than the high RS group, while the tumor purity was the opposite (P <0.0001) (Figure 6A-6D), which demonstrated the reliability of our model and the strong correlation between the model and the patients' immune infiltrate.
The TME can be divided into four subtypes: immune-enriched, non-fibrotic (IE); immune-enriched, fibrotic (IE/F); fibrotic (F); and immune-depleted (D).We calculated risk scores for melanoma patients in GSE22153 and TCGA using the prognostic model.It was found that IE and IE/F groups have lower risk scores compared to the F and D groups (Figure 6E, 6F).CIBERSORT analysis was performed for immune infiltration of patients from the high-and low RS group.We found that patients in the low RS group had high expression in plasma cells, T cell CD8 + , T cell CD4 + memory activated, and NK cells activated, while Macrophages M0 and Macrophages M2 were strongly expressed in the high RS group (Figure 7A).Moreover, we also analyzed the correlation between the 6 biomarkers and immune cells, and found that KIR2DL4 was related to T cell CD8 + , T cell CD4 + memory activated, and NK cells activated (Figure 7B).
According to the results of ESTIMATE, TME and CIBERSORT, the difference of prognosis between the above two groups might be closely associated with T cells.Thus, we used the TIDE to evaluate the therapeutic efficaciousness of immune checkpoint suppression in the two groups.The TIDE score of patients in the high RS group had a better score than those in low one (Figure 8A).The scores of the M2 subtype of tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs) and cancer-associated fibroblasts (CAFs) in the high RS group were higher than the low one (Figure 8B-8D).All results showed that the patients in low RS group had greater prognosis under immune checkpoint inhibition therapy.

Experimental validation
In this study, samples of GSE65904 and GSE19234 were graded according to the prognostic model, and the scoring results were sorted from lowest to highest, and based on the median, the people participated in our research were split into two groups, the high RS group and the low RS group.The survival analysis of the two groups based on survival information suggested that people in the low RS group have significantly better outcomes than people in another group (P < 0.01) (Figure 8E, 8G).The ROC time-dependent curve was applied to verify the prognostic model and all AUCs were greater than or equal to 0.65 (Figure 8F, 8H).These results all validated the effectiveness of our model.

DISCUSSION
Melanoma, as a highly malignant tumor derived from melanocytes, accounts for about 3% of all tumors, mainly occurring in the skin, and a few in the mucosa and viscera [33].Among the most common malignant tumors, the incidence of cutaneous malignant melanoma is the third (6.8%~20%) [34], and during these years, the morbidity and mortality have been still increasing.Compared with other solid tumors, malignant melanoma makes the clinical prognosis worse, with a median survival time of about 6 months and a 1year overall survival rate (OS) of 25% [35].Immune checkpoint inhibitor (ICI) has been presented to significantly boost overall survival rate in persons with advanced-stage melanoma.Ipilimumab (the anti-CTLA-4 antibodies) and nivolumab (the anti-PD-1 antibodies) serve as the routine treatment for advanced melanoma nowadays [36].In this study, six genes (ADCYAP1R1, GPI, IFITM1, KIR2DL4, LIF, and NTS) that exhibited strong correlation with the prognosis of melanoma patients were identified.These genes were utilized to develop a prognostic model that demonstrates reliable predictive value and accuracy.Additionally, the model has the capability to predict the prognosis of ICI   therapy.Moreover, notable differences were observed in the infiltration of immune cells and activation of multiple signaling pathways between the two groups.
In this study, six genes (ADCYAP1R1, GPI, IFITM1, KIR2DL4, LIF, and NTS) were identified as being closely associated with prognosis among a pool of 1793 melanoma-related genes.This selection was achieved through the implementation of Cox univariate regression analysis, Lasso analysis, and Cox multivariate regression analysis, ultimately leading to the establishment of a prognostic model.Subsequently, individuals diagnosed with melanoma were categorized into either a high-risk score (RS) group or a low RS group based on this model, and its efficacy was subsequently validated.Our findings indicate that ADCYAP1R1, GPI, and NTS may contribute to a poor prognosis in melanoma patients, while IFITM1, KIR2DL4, and LIF are more likely to be associated with a better prognosis in individuals with melanoma.These genes have been manifested to be involved in the tumorigenesis in precedent studies.ADCYAP1R1 (ADCYAP Receptor Type I), a Protein Coding gene, involves in related pathways including synthesis and secretion of Aldosterone as well as Signaling mediated by G-protein coupled receptor (GPCR).GO annotations associated with this gene include activities of G protein-coupled receptor and trans membrane signaling receptor [37].The key role of LGR4, a member of the GPCR family, has been demonstrated in tumor immunoregulation, and tumor immunotherapy strategies targeting LGR4 have been proposed and validated.Inhibition of Rspo-Lgr4 switches the polarization of macrophage in order to facilitate checkpoint blockade therapy [38].GPI can encode the protein that is a member of the glucose phosphate isomerase protein family which has been determined as a moonlighting protein, because it is capable to execute distinct functions mechanistically [39].The encoded protein, plays a role as a neurotrophic factor extracellularly, which can promote survival of skeletal motor and sensory neurons, and also as a lymphokine to induce secretion of immunoglobulin [40].According to the auxiliary function of tumor-secreted cytokine and angiogenic factor, it is also as an autocrine motility factor [41]. Neuromedin N and neurotensin have a common precursor, which is encoded by NTS [42].Neurotensin is a secreted tripeptide extensively distributed in the central nervous system and it may be involved in maintenance of intestinal structure and function, and regulation of fat metabolism [43].Diseases associated with IFITM1 include Influenza and Dengue Virus [44].Its related pathways are Interferon gamma signaling [45] and Innate Immune System [46].GO annotations connected with this gene include obsolete signal transducer activity, downstream of receptor [47].Killer cell immunoglobulin-like receptors (KIRs) are expressed by NK cells and subsets of T cells as transmembrane glycoproteins [48].KIR proteins are deemed to play a part in regulation of the immune response because the ligands belonging to subsets of HLA class molecules [49].Leukemia and Myeloid Leukemia are associated with LIF (Interleukin 6 Family Cytokine), which can code related protein [50].Among its related pathways is Interleukin-6 family signaling [51].Therefore, the above studies suggested the rationality and feasibility of the screened 6 genes in the prediction of melanoma incidence and prognosis.
The analysis of differential expression on the genes of patients in the two groups was completed, resulting in the identification of 14 genes highly expressed in the high RS group and 249 genes in the other group.The results of pathway enrichment analysis demonstrated that these genes were mainly associated with several pathways including lymphocyte activation regulation, activation of T cell, cytotoxicity of NK cell, and differentiation of Th1 and Th2 cell.In contrast to T cells, which exhibit specificity towards a singular aberrant molecule on cancerous cells and initiate a targeted assault, natural killer (NK) cells have demonstrated their versatility in the initial defense against cancer.NK cells not only serve as crucial effector cells in the innate immune response, but also function as regulators of diverse immune cell populations.Joy Hsu et al. [52] found that the level of PD-L1-positive NK cells was specifically connected with the outcome of patients with cancer after analyzing the NK cells in human and mice (tumor cells were PD-L1-negative).Under the action of immune checkpoint inhibitors, PD-L1positive NK cells can not only eliminate tumor cells immediately, but also secrete cytokines to control tumor growth.Combined with NK cell activating factor on the basis of PD-L1 antibody can also significantly improve the therapeutic effect.Therefore, NK cells still have more potential to be explored in tumor therapy.Previous studies have shown that in tumor microenvironment (TME), immune infiltration is vital in tumor genesis and progression, and affects the clinical prognosis of tumor patients [53,54].In this study, we found that the melanoma prognostic model constructed using these 6 genes was highly correlated with the patient's immune infiltration.To explore the correlation with circumstance on immune cell infiltration deeply, we continued to use CIBERSORT analysis to compare patients in the two groups, and found that different groups expressed different immune cell subtypes.People in the low RS group showed high expression in plasma cells, T cell CD8 + , T cell CD4 + memory activated, and NK cells activated, while Macrophages M0 and Macrophages M2 were highly expressed in another group.Ali et al. [55,56] have shown that the imbalance of the proportion of immune cell components has strongly correlation with poor outcome and low survival rate of patients with cancer.The past studies have covered that T cell CD8 + and NK cells perform a vital role in tumor immunity [57].Furthermore, an examination was conducted to investigate the association between the aforementioned 6 biomarkers and immune cells.The results revealed a significant correlation between KIR2DL4 and T cell CD8 + , T cell CD4 + memory activated, as well as NK cells activated.It is suggested that the abundant expression of T cell CD8 + and NK cells might reduce the risk elements related to melanoma and ameliorate the prognosis of sick people.
The ICI therapy is the main treatment in the current therapeutic method of melanoma which can significantly upgrade the overall survival rate of people with advanced melanoma [36].CAFs, MDSCs and TAMs can restrict infiltration of T lymphocyte in tumors [58].In the study, we found that the scores of CAFs, MDSCs and TAMs in the high RS group were better than those in another group, which means the people with high risk have higher T cell exclusion and lower infiltration of cytotoxic T lymphocytes.Although the ICI enhance the survival rate of skin melanoma patients, only some patients could be benefited from the ICI treatment [59].Meanwhile, patients are also faced with heavy economic and psychological burden in the treatment process.The TIDE score could presage the response to ICI therapy [60].A positive correlation is observed between lower TIDE scores and improved prognosis.Our study reveals that the low-risk group exhibits lower TIDE scores, indicating potential benefits from ICI therapy.Consequently, the model may serve as a valuable tool for determining the suitability of ICI treatment.
Our research inevitably presents opportunities for further improvement.Firstly, the study entails a bioinformatics analysis utilizing publicly available databases, thereby limiting the authenticity of the molecular mechanism analysis results due to the absence of in vivo or in vitro experiments.Secondly, the sample data obtained from the public database is subject to restrictions, potentially introducing random errors.Moreover, our study only provides some new potential research targets for the prognosis and treatment of melanoma.But the deeper molecular mechanism still needs to be further explored.

CONCLUSIONS
In our study, an effective prognostic model for melanoma was established.All bioinformatics results demonstrated the potential of the six key genes as prognostic markers of melanoma patients and the strong correlation of the model with immune infiltration.This study has provided more information for the pathogenesis and clinical treatment of melanoma.

Figure 2 .
Figure 2. Risk score and ROC curve of melanoma.(A) Risk score of melanoma patients distributed in ascending order.(B) Survival time and status of melanoma patients in order of increasing risk score.The red dots represent the surviving patients and the blue dots represent dead.(C) The heatmap shows the expression of these six biomarkers in melanoma in order of increasing risk scores.(D) The ROC curve for 1, 2, 5-year survival prediction with AUC value.

Figure 4 .
Figure 4. Identification of DEmiRNAs and ceRNA network.(A) Volcano plot of miRNAs between high RS and low RS groups.(B) Heatmap plot of miRNAs between high RS and low RS groups.(C) Differential expression of miRNAs-mRNAs network in melanoma.

Figure 5 .
Figure 5. Somatic mutation analysis.(A) Boxplot of VAF in high RS group.(B) Boxplot of VAF in low RS group.(C) Forestplot of mutant genes in different groups.

Figure 6 .
Figure 6.ESTIMATE analysis.(A) Immune score of high RS samples and low RS samples by ESTIMATE.(B) Stromal score of high RS samples and low RS samples by ESTIMATE.(C) ESTIMATE score of high RS samples and low RS samples by ESTIMATE.(D) Tumor purity of high RS samples and low RS samples by ESTIMATE.(E) Risk scores for four TME subtypes in samples of GSE22153.(F) Risk scores for four TME subtypes in melanoma samples of TCGA.

Figure 7 .
Figure 7. CIBERSORT analysis.(A) The heatmap of immune cell infiltration in high RS samples and low RS samples.Red represents low RS samples, and blue represents high-risk samples.(B) Correlation analysis between different immune cells and biomarkers.

Figure 8 .
Figure 8. TIDE analysis and validation.(A) TIDE score of high RS samples and low RS samples.(B) CAF score of high RS samples and low RS samples.(C) MDSC score of high RS samples and low RS samples.(D) TAM M2 score of high RS samples and low RS samples.(E) Survival analysis of risk score in GSE65904.(F) The ROC curve for 4, 6, 8, 10-year survival prediction with AUC value.(G) Survival analysis of risk score in GSE19234.(H) The ROC curve for 10, 20, 30-month survival prediction with AUC value.