Identification of an Immune-Related Prognostic Signature Associated With Immune Infiltration in Melanoma

Abstract

Melanoma is the leading cause of cancer-related death among skin tumors, with an increasing incidence worldwide. Few studies have effectively investigated the significance of an immune-related gene (IRG) signature for melanoma prognosis. Here, we constructed an IRGs prognostic signature using bioinformatics methods and evaluated and validated its predictive capability. Then, immune cell infiltration and tumor mutation burden (TMB) landscapes associated with this signature in melanoma were analyzed comprehensively. With the 10-IRG prognostic signature, melanoma patients in the low-risk group showed better survival with distinct features of high immune cell infiltration and TMB. Importantly, melanoma patients in this subgroup were significantly responsive to MAGE-A3 in the validation cohort. This immune-related prognostic signature is thus a reliable tool to predict melanoma prognosis; as the underlying mechanism of this signature is associated with immune infiltration and mutation burden, it might reflect the benefit of immunotherapy to patients.

Introduction

Melanoma is the most fatal form of skin tumors and is malignantly transformed from melanocytes (Jackett and Scolyer, 2019). Despite awareness that melanoma progression is primarily caused by complex interactions between environmental and host risk factors (Schadendorf et al., 2018), there were approximately 96,000 new cases and 9,000 deaths of melanoma worldwide in 2018 (Siegel et al., 2018). It is one of the most aggressive cancers that cause more than 75% of skin cancer-related deaths (Gershenwald and Guy, 2015). Although standard surgical resection has improved the prognosis of melanoma patients at localized stage, a large subset of melanoma patients diagnosed with advanced or metastatic melanoma, remains poorly treated.

Due to the complexity of cancer immune response, the contexture and components of immune infiltrates are considerably heterogeneous among each patient and cancer type (Angell and Galon, 2013). Increasing evidence demonstrates that melanoma initiation, evolution, and metastasis are closely related to the tumor microenvironment (TME). In fact, the density of immune contexture has demonstrated a clear correlation with responses to immunotherapy or patients’ clinical outcomes including melanoma (Mann et al., 2013), ovarian clear cell carcinoma (Tan et al., 2019), and bladder cancer (Song et al., 2019). However, the limitations of immune-checkpoint inhibitors (ICIs) in metastatic melanoma require a deep understanding of the definitive mechanisms within the tumor microenvironment and a more reliable biomarker to predict patient prognosis or response to ICIs.

Currently, genetic studies suggest that melanoma is one of the highest mutated malignancies (Alexandrov et al., 2013), of which two prominent mutational events, BARF and NRAS, are observed in approximately 40–60% (Flaherty et al., 2010) and 15–20% (Mandala et al., 2014) of all melanoma cases, respectively. Clinical implications have been observed between mutational burden and susceptibility to immune treatment in the field of oncology (Rizvi et al., 2015; Hugo et al., 2016).

Melanoma, serves as an immunogenic malignancy, developed on the background of these cellular mechanisms including complex interaction among immunomodulatory molecules and high tumor mutational burden (TMB). Given the circumstance that a higher mutation pattern is not always equal to an immunologically hot phenotype (McGranahan et al., 2016), an improved prognostic signature that simultaneously considers the TME and TMB is urgently needed.

In the present study, we downloaded gene expression profiles from The Cancer Genome Atlas (TCGA) as a training cohort and those from the Gene Expression Omnibus (GEO) database as validation cohorts. Based on an immune profile, we used two computational algorithms to screen characteristic genes for further constructing an immune signature with predictive power. We then highlighted the prognostic value and independent role of the resulting 10-gene immune-related model. We also estimated the TME infiltration and mutation pattern in patients from high- and low-risk groups. As a result, we established a robust prognostic biomarker with significantly different TME infiltration and TMB patterns, which can be a potential tool for immunotherapy prediction.

Materials and Methods

Screening DE IRGs Between Primary and Metastatic Melanoma

Using the list of immune genes downloaded from the ImmPort website (Bhattacharya et al., 2018)¹, we screened their expression data in TCGA SKCM data matrix profile during the development of primary and metastatic melanoma. A total of 435 melanoma samples were included in our analysis and the criterion and data processing could be found in our previous study (Liu et al., 2019). Based on |log₂FC| >1 and false discovery rate (FDR) <0.01, differentially expressed immune-related genes (DE IRGs) were obtained using the edge R package. Subsequently, two feature selection methods, the Least Absolute Shrinkage and Selector Operation (LASSO) (Tibshirani, 1996) and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) (Huang et al., 2014), were applied for in-sample cross-validation and reducing the scope of candidate genes for patients with melanoma. Finally, genes from either the LASSO or SVM-RFE algorithms were subsumed in further analysis.

Definition of an Immune-Related Prognostic Model

The feature genes were analyzed using univariate and multivariate Cox proportional hazards regression analysis to establish a prognostic predictive model. A risk score formula was constructed as described in our previous study (Liu et al., 2019). To validate the prognostic capability of the model in TCGA set, melanoma patients were divided into high- or low-risk groups determined by the median score as the threshold. Time-dependent receiver operating characteristic (ROC) curve analyses were utilized to evaluate the accuracy and efficiency of the prognostic model using the “survival ROC” package of R. P < 0.05 was considered statistically significant. The Kaplan-Meier method was used to compare significant differences in overall survival (OS) between different subgroups. Correlations between the risk score and clinical features of melanoma patients were analyzed using the Chi-square test.

Development of a Predictive Nomogram

Univariate Cox regression analysis was performed to identify the independent prognostic factors. The clinical characteristics reaching the statistical difference with P < 0.05 in univariate analysis could be selected for further nomogram construction. Predictive nomograms were constructed using logistic and Cox regression analysis, respectively. Moreover, calibration curves for 5- and 10-year OS were used to assess the prognostic accuracy of the nomogram. The observed and predicted outcomes of the nomogram were presented in the calibrate curve and the 45° line represented the ideal prediction. The nomograms and calibration curves were constructed using the “rms” package in R.

Risk-Related Genes Analysis

To explore the biological implications of the immune-related prognostic signature, the melanoma patients were divided into low- and high-risk groups to identify the differentially expressed genes (DEGs) associated with the risk patterns mentioned above. DEGs among these groups were determined using the R package “limma,” with |log₂FC| >1 and FDR <0.01 regarded as significance criteria.

Gene ontology (GO) enrichment analysis within risk groups was performed to gain insight into the biological functions of prognostic immune-related genes using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Huang et al., 2009), with a threshold of P < 0.05. Furthermore, significantly enriched GO terms in biological process (BP) were visualized using R packages. Gene set variation analysis (GSVA) (Hänzelmann et al., 2013) was conducted between high- and low-risk groups using the GSVA_1.30.0 package of R.

Estimation of Tumor Infiltrating Immune Cells

To investigate the immune infiltration landscape between different subgroups, single-sample gene set enrichment analysis (ssGSEA) was used to quantify the infiltration levels of 24 immune cell types, according to the expression levels of the immune cell-specific marker genes described by Bindea et al. (2013). The ssGSEA ranked the genes by their absolute expression in each sample and grouped the different immune cell infiltration patterns based on the R package “gsva.”

Similarly, cell type identification by estimating the relative subsets of RNA transcripts (CIBERSORT)², a deconvolution algorithm, was developed to determine the relative fraction of 22 immune cell types in melanoma tissues (Newman et al., 2015). The transcriptional profiles obtained from TCGA database were prepared in accordance with the accepted format of CIBERSORT and LM22, a gene signature matrix that defines 22 immune cell subtypes, and was used as the signature gene file. CIBERSORT was run with 1,000 permutations and a threshold <0.05 as recommended.

Based on the gene expression signature, “Estimation of STromal and Immune cells in Malignant Tumors using Expression data” (ESTIMATE) was used to assess the estimate score, stromal score, and immune score using the “estimate” package for tumor samples (Yoshihara et al., 2013)³.

Profiles of TMB and Correlation Analysis

The mutation data for melanoma patients were downloaded from TCGA data portal⁴. For each patient, MutSigCV_v1.41⁵ was used to identify the significant mutated genes (P < 0.05) across the two classes currently identified with the risk score (Lawrence et al., 2013). The mutation landscape oncoprint was drawn using R package “ComplexHeatmap.”

Validation Datasets

The robustness of the prognostic model was validated in GSE19234, GSE22153, and GSE35640 based on the GPL570, GPL6102, and GPL570 platforms from the GEO database (Clough and Barrett, 2016)⁶, respectively.

Availability of Data and Materials

All data used in this study were obtained from the TCGA database: https://portal.gdc.cancer.gov/and GEO database: https://www.ncbi.nlm.nih.gov/geo/.

Statistical Analysis

All statistical analyses were executed in R version 3.5.2. P < 0.05 were considered statistically significant unless otherwise mentioned.

Results

Screening Candidate DE IRGs

In total, 288 DE IRGs were screened out with the cutoff criteria of |log₂FC| >1 and FDR <0.01, including 221 upregulated and 67 downregulated genes. For further validation and selection of DE IRGs with significantly characteristic values to classify primary melanoma (PM) and metastatic melanoma (MM), the LASSO algorithm was used to identify a set of 190 DEGs (Figures 1A,B) and the SVM-RFE algorithm was used to select a set of 55 DEGs (Figures 1C,D). Finally, after combing the LASSO and SVM-RFE algorithms, the 207 most representative DEGs were identified, with 38 genes selected simultaneously by the two algorithms (Figure 1E).

FIGURE 1

Construction of a Prognostic Immune Gene Signature

In this study, each DE IRG was first submitted for univariate Cox proportional hazards regression with the criteria of a P < 0.05 (Supplementary Figure S1A). Then, multivariate Cox regression analysis was performed to develop the prognostic model and the Concordance index was 0.7, which indicated the high predictive accuracy of the signature for survival (P < 0.05; Table 1 and Supplementary Figure S1B). The prognostic risk score was calculated as follows: (0.182 ^∗ S100A12) + (0.390 ^∗ CCRL2) + (−0.587 ^∗ CD86) + (−0.2 90 ^∗ IL21R) + (0.199 ^∗ CCR4) + (0.375 ^∗ FCGR3A) + (−0.350 ^∗ KLRD1) + (−0.147 ^∗ IL18RAP) + (0.460 ^∗ IL2RB) + (−0.2 11 ^∗ CCL8).

Melanoma patients were assigned to the high- and low-risk groups based on the risk score model. Distribution of risk score, patients’ survival status, and gene expression profiles associated with the risk score are shown in Figure 2A. Kaplan-Meier analysis demonstrated that melanoma patients with a high-risk score showed dramatically worse prognosis than those with a low-risk score (Figure 2B), and the area under the curve values of the time-dependent ROC curve were 0.731, 0.774, and 0.76 for 3-, 5-, and 10-year survival, indicating the high specificity and sensitivity of the prognostic signature (Figure 2C).

FIGURE 2

Stratification Analyses of the 10-Gene Prognostic Signature

Considering the potential impacts of clinical characteristics on the risk score of the prognostic model, a Chi-square test was performed. The risk score exhibited a higher level in event and the American Joint Committee on Cancer-Tumor (AJCC-T) subgroups; there were no significant differences in the risk score, AJCC-Nodes (AJCC-N), AJCC-Metastasis (AJCC-M), Stage, Clark, Primary Tumor, and the BRAF and NRAS subgroups (Supplementary Figure S2A). To assess whether the prognostic classifier was an independent indicator in distinct subgroups, we checked the effect of BRAF and NRAS mutation or wild-type on the prognostic ability for melanoma in TCGA cohort. Kaplan-Meier analysis revealed that patients in the high-risk group were associated with a higher mortality risk in the BRAF^wt and NRAS^wt or BRAF^mut and NRAS^mut subgroups (Supplementary Figures S2B–E). Additionally, patients in stages I and II or stages III and IV were divided into two groups with the median value. We found that patients with high risk were classified into the group with a worse prognosis, compared to patients with a low risk (Supplementary Figures S2F,G). These results demonstrated the robust and predictive power of this prognostic model, in which patients high-risk scores had a shorter overall survival than those with low risk scores in each stratum.

Construction of a Nomogram Model

To develop a clinically relevant quantitative approach for predicting the survival probability of a patient with melanoma, we constructed predictive nomograms. Based on univariate analysis (Figure 3A), we generated two nomograms to predict the death odds of patients using logistic regression and survival rate with Cox regression analysis (Figures 3B,C). The calibration plots for the 5- and 10-year survival showed an optimal agreement between the nomogram-predicted and observed OS (Figure 3D).

FIGURE 3

Functional Traits of the Prognostic Signature

To explore the potential cause of the prognostic signature, we divided the melanoma patients from TCGA database into high- and low-risk groups, based on the median risk score. After edgeR filtering (|log₂FC| > 1 and FDR < 0.05), we screened out 1,251 DEGs, including 552 upregulated and 699 downregulated genes (Figure 4A). Of these DEGs, 26 genes were immune-related and are highlighted in Figure 4A. GO enrichment analysis revealed that upregulated genes were significantly enriched in multiple pathways including T cell activation, regulation of T cell activation, and regulation of lymphocytes (P < 0.05; Figure 4B). Moreover, downregulated genes were significantly enriched in epidermis development, skin development, and keratinization (P < 0.05; Figure 4C). GSVA showed that patients with low risk scores exhibited increased expression of proteins associated with the interferon gamma response, allograft rejection, and interferon alpha response (Figure 4D). These findings indicate differences in the immune-related genes and signaling pathways between high- and low-risk groups, which might partly explain the reason for the significant difference in prognosis between the subgroups.

FIGURE 4

The Risk Score Was Associated With Immune Cell Infiltration

Immune cell infiltration status was assessed by applying the ssGSEA approach to the melanoma transcriptomes. Twenty-four immune-related terms were incorporated to deconvolve the abundance of diverse immune cell types in melanoma. The whole cohort was clustered into two clusters in terms of immune infiltration by applying the IRGs signature (Figure 5A) and the relative immune score in ssGSEA is shown in Supplementary Figure S3A. Subsequently, a TME cells network was constructed with three main clusters depicting a comprehensive landscape of tumor-immune cell interactions and cell lineage, and their effects on the OS of melanoma patients (Supplementary Figure S3B and Figure 5B). ssGSEA was used to assess the relative proportions of the 24 immune cells in each CC sample. The abundance of aDC, B cells, CD8 T cells, cytotoxic cells, DC, macrophages, NK CD56^dim cells, pDC, T cells, T helper cells, Tcm, TFH, Th1, and Treg cells was low in the 10-IRG signature high-risk group and associated with better OS. The abundance of eosinophils and mast cells was high in the 10-IRG signature high-risk group and associated with poor OS.

FIGURE 5

The immune infiltration in melanoma tissues between the high- and low-risk groups was then investigated using the CIBERSORT algorithm. The proportion of 22 immune cells in each subgroup is shown in a bar plot (Figure 6A). The results revealed that plasma cells, CD8 T cells, CD4 memory activated T cells, follicular helper T cells, and M1 macrophages were negatively correlated with the risk score (Figures 6B,C) and that M0 macrophages, M2 macrophages, activated DCs, and neutrophils were positively correlated with the risk score (Figures 6B,C). Figure 6D indicated the poor correlation coefficient between 22 immune cells. The population with negative relation included M0 macrophages and CD8 T cells (r = −0.61), and CD4 memory resting T cells and CD8 T cells (r = −0.57). Further, neutrophils and activated mast cells had a positive relation (r = 0.71). After the ESTIMATE algorithm was processed, a higher Estimate score was found in the low-risk group. Similarly, the fraction of immune and stromal cells was associated with the low-risk group (Figure 6E).

FIGURE 6

Different Tumor Mutation Burden (TMB) Patterns Between Two Risk Groups

We defined and calculated the TMB variable between low- and high-risk groups. We also assessed the correlation between risk score and mutated genes and the mutant rate in these 47 mutants, which were distributed in more than 10% of melanoma samples, and were significantly associated with risk scores at P < 0.05 (Figure 7A). The mutational landscape indicated that mutation events occurred more frequently in the low-risk group than in the high group (Figures 7B,C). Moreover, we further analyzed the survival significance of TMB in melanoma and found that lower TMB levels were associated with worse overall survival outcomes (Figure 7C).

FIGURE 7

Validation of Signature

To substantiate the stability of the prognostic signature, three external validation cohorts were analyzed. For the external validation cohort 1 and 2, GSE19234 and GSE22155 contained 44 and 57 melanoma patients, respectively. Consistent with previous results, the low risk group had higher levels of immune cell infiltration (Supplementary Figures S4A,S5A). The relative immune score is shown in Supplementary Figures S4B,S5B. As expected, patients in the high-risk group had a significantly increased mortality risk compared with those in the low-risk group (Supplementary Figures S4C,S5C). For the external validation cohort 3, GSE35640 contained 65 patients with MAGE-A3 antigen-specific cancer immunotherapy. Similar analysis showed that the low-risk group had higher levels of immune cell infiltration (Supplementary Figure S6). Furthermore, the prognostic model demonstrated the potential to predict the effect of MAGE-A3 immunotherapy in melanoma patients (Supplementary Figure S6).

Discussion

The clinical heterogeneity of melanoma suggested that biologically relevant differences might exist within subtypes. The purpose of the study was to identify an immune gene expression signature to predict the prognosis in patients between primary and metastatic melanoma. A flowchart of the analysis procedure for this study was shown in Supplementary Figure S7. According to the DE IRGs obtained in our study, we constructed a 10-IRG prognostic signature independent of previously known clinicopathological factors such as the mutant status. The enriched terms for biological process and immune infiltration landscape revealed that the TME might contribute to tumor progression and a poorer prognosis in melanoma, which was validated in two external datasets. Subsequently, our data demonstrated that most tumors could be categorized as either high TMB with a better prognosis or low TMB with a worse prognosis. This classification of melanoma identified in our study demonstrates that patients with a low risk score had increased immune cell infiltration and TMB, and were more likely to respond to ICIs.

We constructed and validated an immune-related risk signature for melanoma using TCGA and GEO datasets. The signature was composed of 10 DE IRGs with prognostic capability. In this prognostic model, five DE IRGs (S100A12, CCRL2, CCR4, FCGR3A, and IL2RB) were used as risk factors with positive coefficients, whereas the other five genes (CD86, IL21R, KLRD1, IL18RAP, and CCL8) were protective factors with negative coefficients. Furthermore, there were significant differences in survival curves between patients with high and low risk scores. The signature exhibited a firm predicting capability in the training and validation datasets and the high prognostic categorization performance of the immune signature was assuredly due to our stratified analysis strategy. Furthermore, two nomograms were built based on univariate Cox regression coefficients of Risk, Clark status, Stage, AJCC-N, and AJCC-T to evaluate the prognosis predictions of the immune signature. A satisfactory agreement between the observed values and the predicted values for the 5- and 10-years OS, was observed in the calibration curves. These results indicated that this model might be an effective tool for outcome prediction of individualized melanoma patients.

Moreover, we reclassified the microarray data into DEGs according to the median risk score. Functional enrichment analysis indicated that risk-related DEGs were primarily involved in a multitude of immune pathways. We speculated that locoregional immune status might have the potential to improve melanoma prognosis classification. The TME consists of cellular components (immune cells, etc.) as well as non-cellular components (cytokines, etc.). Notably, the complex interplay between tumor cells and their surrounding microenvironment plays a pivotal role during tumor development and has significant effects on the OS and immunotherapeutic efficacy in tumors (Mann et al., 2013; Song et al., 2019; Tan et al., 2019). Here, we investigated the immune landscape between patients in the high-risk and low-risk groups using two bioinformatic methods to infer specific immune cell infiltration. In our analysis, melanoma patients were clustered into two main clusters (high or low immune infiltration). Patients with better prognosis were found in the high infiltration cluster. Furthermore, infiltration of CD103 + CD8 + T cells has been associated with longer survival in patients with melanoma tumors (Edwards et al., 2018). Among these immune cells, the abundance of CD8 T cells was high in the 10-IRG signature low-risk group and was associated with better OS based on the results of the ssGSEA and CIBERSORT algorithms. Based on the results of the ssGSEA, the abundance of immune cells with anti-tumor effects including CD8 T cells, TH1 cells and B cells, was high in low-risk group, though the distribution in high-risk group showed less immune cells, named “cold tumor” (Galon and Bruni, 2018). Subsequently, CIBERSORT was used to assess the relative proportions of 22 immune cells in each melanoma sample. The abundance of immune cells with anti-tumor effects including CD8 T cells, CD4 T cells and M1 macrophages, was high in the low-risk group. However, in high-risk group, the abundance of M2 macrophages required for cancer progression and unactivated M0 macrophages were high, indicating that there was a tendency to transform into M2 macrophages with anti-inflammatory and tumor promoting effects (Mills et al., 2016). And previous studies have been reported that patients with more anti-tumor immune cells infiltration, named “hot tumor,” would have better survival prognosis (Galon and Bruni, 2018). Furthermore, using the ESTIMATE method to evaluate the main non-tumor components in the microenvironment, we found that the low-risk group was profoundly associated with either a higher immune score or a higher stromal score. These results may partially explain the predictive value of this signature.

Due to the prevalence of somatic mutations in the melanoma genome, further mutation analysis was performed to explore the possible mechanisms underlying the signature’s prognostic value. In our model, the risk score was in contrast to the TMB pattern to determine the prognosis of melanoma patients, suggesting that the poor prognosis of the high-risk group may be due to fewer mutant genes in this group. Recently, some studies have indicated that increased tumor mutation loads were associated with survival benefit from both anti-CTLA-4 and anti-PD-1 therapy in multiple malignancies such as melanoma (Hugo et al., 2016), lung cancer (Rizvi et al., 2015), and esophagogastric cancer (Greally et al., 2019). Consistent with these studies, our results showed that melanoma patients in the low-risk group with higher immune infiltration and tumor mutation load might respond well to immunotherapy, though larger studies will be necessary to confirm this finding.

Previous studies have explored whether immune-related biomarkers could be indicators for prognosis or response to immune therapy of melanoma (Cursons et al., 2019; Nie et al., 2019; Poźniak et al., 2019; Yang et al., 2020). Compared with these studies, we used a combination strategy to control the robustness of the predictions of DE IRGs. Based on stratified analysis, we also determined the predictability of the nomogram description model, which fully demonstrated the predictability of the model identified in this study. Poźniak et al. (2019) and Yang et al. (2020) evaluated the abundance of infiltrating immune cells between prognostic immune subgroups. However, our study systematically evaluated the landscape of immune cells in melanoma samples using three bioinformatics tools for cross validation. Notably, our prognostic model demonstrated the compatibility of immune signature with TMB in a multi-omics manner, which should be highlighted during the management of immunotherapeutic combinations. Ultimately, we used this model on its immune profile to evaluate the response to ICIs in melanoma and observed positive results in this study.

Despite these promising results, our study has some limitations. The melanoma samples included in our studies were obtained from available public data; the clinical utility of the model needs to be confirmed in a large-scale of melanoma patients. Though the prognostic value of the immune-related signature was validated well in two GEO cohorts, supplemental basic experiments are still warranted to uncover the biological mechanisms of 10 DE IRGs in the promotion of tumor development. Furthermore, due to limited data, the relevance of the signature to ICI response is not fully understood and requires further research.

Conclusion

In the current study, we performed a comprehensive evaluation of the prognostic signature generated and validated in our study, which might be a clinically promising tool to classify melanoma patients into subgroups with distinct outcomes, immunophenotypes, mutation patterns, and even their distinct responses to immune therapy. It provides new implications regarding the melanoma immune microenvironment, TMB, and immune-related therapy.

References

• 1

  AlexandrovL. B.Nik-ZainalS.WedgeD. C.AparicioS. A.BehjatiS.BiankinA. V.et al (2013). Signatures of mutational processes in human cancer.Nature500415–421.

  • Google Scholar
  • View reference in article

• 2

  AngellH.GalonJ. (2013). From the immune contexture to the immunoscore: the role of prognostic and predictive immune markers in cancer.Curr. Opin. Immunol.25261–267. 10.1016/j.coi.2013.03.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BhattacharyaS.DunnP.ThomasC. G.SmithB.SchaeferH.ChenJ.et al (2018). ImmPort, toward repurposing of open access immunological assay data for translational and clinical research.Sci. Data5:180015.

  • Google Scholar
  • View reference in article

• 4

  BindeaG.MlecnikB.TosoliniM.KirilovskyA.WaldnerM.ObenaufA. C.et al (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer.Immunity39782–795. 10.1016/j.immuni.2013.10.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CloughE.BarrettT. (2016). The gene expression omnibus database.Methods Mol. Biol.141893–110. 10.1007/978-1-4939-3578-9_5

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  CursonsJ.Souza-Fonseca-GuimaraesF.ForoutanM.AndersonA.HollandeF.Hediyeh-ZadehS.et al (2019). A gene signature predicting natural killer cell infiltration and improved survival in melanoma patients.Cancer Immunol. Res.71162–1174. 10.1158/2326-6066.cir-18-0500

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  EdwardsJ.WilmottJ. S.MadoreJ.GideT. N.QuekC.TaskerA.et al (2018). CD103+ tumor-resident CD8+ T cells are associated with improved survival in immunotherapy-naïve melanoma patients and expand significantly during anti-PD-1 treatment.Clin. Cancer Res.243036–3045. 10.1158/1078-0432.ccr-17-2257

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  FlahertyK. T.PuzanovI.KimK. B.RibasA.McArthurG. A.SosmanJ. A.et al (2010). Inhibition of mutated, activated BRAF in metastatic melanoma.N. Engl. J. Med.363809–819.

  • Google Scholar
  • View reference in article

• 9

  GalonJ.BruniD. (2018). Approaches to treat immune hot, altered and cold tumours with combination immunotherapies.Nat. Rev. Drug Discov.18197–218. 10.1038/s41573-018-0007-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  GershenwaldJ. E.GuyG. P.Jr. (2015). Stemming the rising incidence of melanoma: calling prevention to action.J. Natl. Cancer Inst.108:djv381.

  • Google Scholar
  • View reference in article

• 11

  GreallyM.ChouJ. F.ChatilaW. K.MargolisM.CapanuM.HechtmanJ. F.et al (2019). Clinical and molecular predictors of response to immune checkpoint inhibitors in patients with advanced esophagogastric cancer.Clin. Cancer Res.256160–6169. 10.1158/1078-0432.ccr-18-3603

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  HänzelmannS.CasteloR.GuinneyJ. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data.BMC Bioinform.14:7. 10.1186/1471-2105-14-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  HuangD. W.ShermanB. T.LempickiR. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.Nat. Protoc.444–57. 10.1038/nprot.2008.211

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  HuangM. L.HungY. H.LeeW. M.LiR. K.JiangB. R. (2014). SVM-RFE based feature selection and Taguchi parameters optimization for multiclass SVM classifier.Sci. World J.2014:795624.

  • Google Scholar
  • View reference in article

• 15

  HugoW.ZaretskyJ. M.SunL.SongC.MorenoB. H.Hu-LieskovanS.et al (2016). Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma.Cell16535–44. 10.1016/j.cell.2016.02.065

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  JackettL. A.ScolyerR. A. (2019). A review of key biological and molecular events underpinning transformation of melanocytes to primary and metastatic melanoma.Cancers11:E2041.

  • Google Scholar
  • View reference in article

• 17

  LawrenceM. S.StojanovP.PolakP.KryukovG. V.CibulskisK.SivachenkoA.et al (2013). Mutational heterogeneity in cancer and the search for new cancer-associated genes.Nature499214–218.

  • Google Scholar
  • View reference in article

• 18

  LiuN.LiuZ.LiuX.ChenH. (2019). Comprehensive analysis of a competing endogenous RNA network identifies seven-lncRNA signature as a prognostic biomarker for melanoma.Front. Oncol.9:935. 10.3389/fonc.2019.00935

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  LiuN.LiuZ.LiuX.DuanX.HuangY.JinZ.et al (2020). Identification of an immune-related prognostic signature associated with immune infiltration in melanoma.Research Square [Preprint], 10.21203/rs.3.rs-21119/v1

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  MandalaM.MerelliB.MassiD. (2014). Nras in melanoma: targeting the undruggable target.Crit. Rev. Oncol. Hematol.92107–122. 10.1016/j.critrevonc.2014.05.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  MannG. J.PupoG. M.CampainA. E.CarterC. D.SchrammS. J.PianovaS.et al (2013). BRAF mutation, NRAS mutation, and the absence of an immune-related expressed gene profile predict poor outcome in patients with stage III melanoma.J. Invest. Dermatol.133509–517. 10.1038/jid.2012.283

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  McGranahanN.FurnessA. J.RosenthalR.RamskovS.LyngaaR.SainiS. K.et al (2016). Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade.Science3511463–1469.

  • Google Scholar
  • View reference in article

• 23

  MillsC. D.LenzL. L.HarrisR. A. (2016). A breakthrough: macrophage-directed cancer immunotherapy.Cancer Res.76513–516. 10.1158/0008-5472.can-15-1737

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  NewmanA. M.LiuC. L.GreenM. R.GentlesA. J.FengW.XuY.et al (2015). Robust enumeration of cell subsets from tissue expression profiles.Nat. Methods12453–457. 10.1038/nmeth.3337

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  NieR. C.YuanS. Q.WangY.ChenY. B.CaiY. Y.ChenS.et al (2019). Robust immunoscore model to predict the response to anti-PD1 therapy in melanoma.Aging1111576–11590. 10.18632/aging.102556

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  PoźniakJ.NsengimanaJ.LayeJ. P.O’SheaS. J.DiazJ. M. S.DroopA. P.et al (2019). Genetic and environmental determinants of immune response to cutaneous melanoma.Cancer Res.792684–2696. 10.1158/0008-5472.can-18-2864

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  RizviN. A.HellmannM. D.SnyderA.KvistborgP.MakarovV.HavelJ. J.et al (2015). Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer.Science348124–128.

  • Google Scholar
  • View reference in article

• 28

  SchadendorfD.van AkkooiA. C. J.BerkingC.GriewankK. G.GutzmerR.HauschildA.et al (2018). Melanoma.Lancet392971–984.

  • Google Scholar
  • View reference in article

• 29

  SiegelR. L.MillerK. D.JemalA. (2018). Cancer statistics.CA Cancer J. Clin.687–30.

  • Google Scholar
  • View reference in article

• 30

  SongB. N.KimS. K.MunJ. Y.ChoiY. D.LeemS. H.ChuI. S. (2019). Identification of an immunotherapy-responsive molecular subtype of bladder cancer.eBio Med.50238–245. 10.1016/j.ebiom.2019.10.058

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  TanT. Z.YeJ.YeeC. V.LimD.NgoiN. Y. L.TanD. S. P.et al (2019). Analysis of gene expression signatures identifies prognostic and functionally distinct ovarian clear cell carcinoma subtypes.eBio Med.50203–210. 10.1016/j.ebiom.2019.11.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  TibshiraniR. (1996). Regression shrinkage and selection via the LASSO.J. R. Stat. Soc. Ser. B58267–288. 10.1111/j.2517-6161.1996.tb02080.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  YangS.LiuT.NanH.WanågY.ChenH.ZhangX.et al (2020). Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of cutaneous melanoma.J. Cell. Physiol.2351025–1035. 10.1002/jcp.29018

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring tumour purity and stromal and immune cell admixture from expression data.Nat. Commun.41–11.

  • Google Scholar
  • View reference in article