Transcriptome (n = 472) and clinical feature data (n = 477) were downloaded from the UCSC dataset.¹ The RNA-seq data included 471 melanoma samples and 1 normal skin sample. The clinical information of melanoma samples included survival events, tumor staging, grade, Breslow depth value, and Clark level value. We selected samples (n = 472) with both RNA-seq data and clinical phenotype data for subsequent analysis. Later, when building the Lasso regression model, we only retained samples (n = 452) containing both RNA-seq data and survival data for analysis. We installed a 2:1 ratio to randomly group these samples as the training set (n = 301) and internal validation set (n = 151), respectively. All the data used for survival analysis comes from an attachment to an article (Liu et al., 2018). The somatic mutation data (MuTect2 Variant Aggregation and Masking) were also downloaded from this website. We built a flow chart for the whole process of the experiment (Figure 1). We converted our downloaded RNA-seq data (FPKM values) into TPM values, which makes the samples more comparable (Wagner et al., 2012). To filter out genes whose expression has not changed much, we then selected the genes in the top 25% (n = 4,939) of the variance for subsequent analysis. For the processing of externally validated data from the GEO database, we download public data from the GEO database, then RMA corrected, standardized, and then logarithmic.

We used those genes from one step to build a co-expression network by the R package “WGCNA” (Langfelder and Horvath, 2008). First, the Pearson correlation coefficient between samples is calculated to remove abnormal samples (Supplementary Figure 1). Then we constructed an adjacency matrix through the use of Pearson’s correlation analysis for all pairs of genes. Finally, by associating the modules with the clinical features, we selected the modules most relevant to the clinical features of greatest concern for subsequent analysis. We defined the cor.geneModuleMembership > 0.8 and cor.geneTraitSignificance > 0.2 as the threshold for screening hub genes in a module.

We uploaded all the genes of the blue module to DAVID’s website² for GO analysis and KEGG analysis. In order to understand the enrichment of signal pathways in the classification of the 5-mRNA signature, the gene set enrichment analysis (GSEA) was carried out with R package “ClusterProfiler.”

The mutation frequency less than 0.05 and synonymous mutations were removed. We used CIBERSORT (Newman et al., 2015) online software to quantify the immune infiltration of BLCA-SKCM based on deconvolution algorithm. The LM-22 gene signature which included 22 types of human immune cell was used as the set of reference gene-expression values to perform immune infiltration analysis. The R package “estimate” (Yoshihara et al., 2013) was used to quantify the stromal and immune cell admixture, stromal cells and immune cells were the main components of normal cells in tumor tissue. The association between tumor stem cell markers (CD20, ABCB5, CD38, etc.) and risk score was calculated using Pearson’s correlation analysis.

We use R package “survival” for univariate and multivariate analysis, R package “forestplot” for forest map, and then R package “rms” for the nomogram. The calibration curve and the decision curve analysis performed by R package “rmda” were used to measure the accuracy of the nomogram. The statistical test method of the two groups of data involved in this paper was the t-test, which used the two-sample t-test and the Welch two-sample t-test, respectively, according to whether the variance was homogeneous. The statistical method used for all survival analyses was the Log-rank test.

The R package of “WGCNA” is used to perform co-expression network analysis. The Pearson correlation analysis was used to detect the outlier samples and we found 8 samples were outlier samples, then we removed them (Supplementary Figure 1). The 4,939 genes with the first 25% variance were put with similar expression patterns into modules by cluster analysis. The power of β = 10 was chosen for the soft-thresholding to develop a scale-free network (scale-free R² = 0.91, Supplementary Figure 2).

We finally got six modules, among which the blue module showed a strong negative correlation with the Breslow depth (r = −0.32, p = 3e-12). It also showed a strong negative correlation with the Clark level (r = −0.22, p = 1e-06) and the T stage (r = −0.25, p = 6e-08, Figure 2A). So, we finally chose the blue module as our focus, and then finally we got 116 hub genes (Figure 2B).

We constructed a prognosis signature based on 5-mRNA through LASSO algorithm (Figures 2C,D). In the training set (HR = 2.69, 95%CI: 1.94–3.71, p < 0.0001) and the external validation set (HR = 1.73, 95%CI: 0.96–3.10, p = 0.047), the prognostic model showed high predictive power. The risk score = −0.08 × SRGN expression value—0.04 × IDO1 expression value—0.04 × RARRES3 expression value—0.03 × CCL4 expression value—0.007 × HLA-DPB1 expression value. We grouped the population according to the median risk score. We found that the risk of death in the training set and the validation set increased significantly with the increase of the risk score (Figures 3A,D, the middle panels), and the genes used to construct the 5-mRNA prognosis signature decreased with the increase of the risk score (Figures 3A,D, the lower panels). The time-dependent ROC curve showed that the AUC of 1-, 3-, 5-, and 7-years was mostly greater than 0.65 (Figures 3B,E). Survival analysis also showed good predictive value of the signature; the high-risk group tended to have a worse prognosis (Figures 3C,F).

Univariate and multivariate Cox regression analyses were performed to test the independence of the predictive power of risk scores. First, the signature risk score, TMN stage, gender, histologic stage, Breslow depth, and Clark level were performed with the univariate analysis (Figure 4A and Table 1), we found that risk score (HR = 5.05, p = 2.91e-11), Breslow depth (HR = 1.03, p = 6.98e-05), Clark level (HR = 1.74, p = 7.19e-07), T stage (HR = 1.33, p = 7.68e-07), N stage (HR = 1.31, p = 7.49e-05), and pathological stage (HR = 1.44, p = 1.42e-05) showed good predictive power of prognosis. Next, we put this factor to perform the multivariate analysis. The results showed that risk score was an independent predictor of melanoma risk (HR = 4.6, p = 8.77e-07, Figure 4B and Table 1). The subgroup analysis’s results exhibit the 5-mRNA signature showed good predictive value among different tumor stages and grades (Supplementary Figure 3). We also used the GEPIA database (Tang et al., 2017)³ to measure the prognostic value of the 5 genes in the signature. The results showed that all the 5 genes were protective factors for the prognosis of melanoma (Supplementary Figure 4).

From the results of the previous step, we constructed a nomogram combining the risk score and important clinical characteristics (Figure 5A). The calibration curve showed that the predictive power of the nomogram at 1, 3, and 5 years was accurate (Figures 5B–D). The DCA indicated the clinical utility of the nomogram was well, the risk score had a better net benefit compared to tumor stage, and the combination of the two had a better net benefit than either of the two (Figure 5E).

We uploaded all the genes of the blue module to DAVID’s website to perform GO and KEGG analysis, and the results of KEGG analysis indicated that the “Staphylococcus aureus infection,” “Antigen processing and presentation,” “Graft-vs.-host disease” were enriched in the blue module (Supplementary Figure 5D and Supplementary Table 1). The GO analysis showed that the blue module was related to “immune response,” “type I interferon signaling pathway,” “inflammatory response,” “T cell receptor signaling pathway,” “extracellular exosome,” “peptide antigen binding,” “MHC class II receptor activity,” and so on (Supplementary Figures 5A–C and Supplementary Table 1). These results indicated that the genes of the blue module were significantly related to the immune process. Then we explored the difference in signaling pathways between the high- and low-risk groups by GSEA under the threshold of p.adjust < 0.05. The high-risk groups were mainly enriched in “MYC targets v2” related pathways, while the low-risk groups were mainly enriched in “ALLOGRAFT REJECTION,” “INTERFERON_ALPHA_RESPONSE,” “INFLAMMATORY_RESPONSE,” and so on (Supplementary Figure 6C and Table 2). Through GO and KEGG analysis of risk scores, we found that risk scores were significantly correlated with T cell activation, regulation of lymphocyte activation, and regulation of T cell activation functions, and were correlated with cell adhesion molecules (CAMs), allograft rejection, antigen processing and presentation and so on pathways (Supplementary Figures 6A,B and Supplementary Table 2).

Through the functional annotation in the previous step, the results exhibit that the low-risk group had a significant relationship with tumor immune-related pathways, so we then analyzed the characteristics of this prognostic model in the immune microenvironment (Figures 6A,B). We quantified 22 major immune cells using CIBERSORT software and found that the main immune infiltration cells in melanoma are CD8 T cells and macrophages (Figures 6A,D and Supplementary Figures 7A–C), the low-risk group had higher CD8 T cell and macrophages M1 infiltration and lower macrophage M0/2 infiltration (Figure 6C). It was also found that these recognized immune checkpoints (PD-L1, PD-L1, CTLA4, etc.) were highly expressed in the low-risk group (Supplementary Figures 7E–G). Our study found that in patients with high TMB or BRAF mutation, the risk score was low (Supplementary Figures 7D,H). By calculating risk score and the correlation of the tumor stem cell marker, results showed that the risk score significantly associated with melanoma cancer stem cell markers (CD20, ABCB5), not only such, risk score also some correlation with other tumor markers, which suggests that tumor stem cell factors contributed to the prognosis of different risk groups for different reasons (Figure 6D).