3.1 Clinical characteristics of the study populations
A total of 211 WT Bf-CM patients’ samples was downloaded from TCGA. As is shown in Table 1, 64.0% of the samples were collected from male and 36% were from female. The median age of this cohort is 62, and the range of age is 25~87. 193 samples were matched with pathological stages according to the American Joint Committee on Cancer (AJCC) Cancer staging manual, the percentages of stage 0, I, II, III, IV were 1.9%, 16.6%, 29.9%, 39.3% and 3.8%, respectively. 25.5% of these samples featured with a Breslow thickness < 2mm, 28.5% within the range of 2mm to 5mm, 25.1% > 5mm. More information of the selected clinical characteristics of are listed in Table 1.
3.2 Construction and validation of the expression-based prognostic signature
In order to generate the training set and test set, 211 WT Bf-CM patients were divided equally into two set at random, and their clinical characteristics are shown in Table 1. By analyzing the transcriptome data of training set with univariate Cox regression, 251 RNA were shown significantly to correspond with the 5-year OS of WT Bf-CM patients (p<0.01). Next, multivariate Cox stepwise regression was performed based on the results of random forests and an expression-based prognostic signature was constructed. This prognostic signature is based on the hazard ratio of 7-RNA and can be formulated as: Risk scores = 0.318*expression of CDC73 + 0.282*expression of RP1-69E11.3 + 0.882*expression of RP11-188D8.1 + 0.968*expression of RP11-116P24.2 + 0.486*expression of TRIB2 + 0.203*expression of VPS13D + 0.199*expression of CELF3. More information about the 7 RNA is listed in Table 2.
The risk scores of each patient from training set were calculated and the median score was used as cut-off for the stratification of high risk group and low risk group. Kaplan–Meier survival analysis and Log-rank test of this two groups showed a significantly poorer prognostics of high risk group compared to low risk group (39.8 months vs. 170.2 months, p<0.001) (Figure.1a). For the prognostic value of the proposed signature, ROC curve of training set was drawn and AUC was calculated to be 0.866 (Figure.1a). To further validate the performance of the proposed signature, patients from test set were stratified similarly based on the median score from training set. Same analysis and tests were also performed with the high- and low-risk groups of test set. Besides the significant correspondence between risk-score and prognosis (59.3 months vs. 90.4 months, p=0.005), the ROC curve (AUC = 0.699) reveals the considerable potential of this signature for predicting the prognosis of WT Bf-CM patients in the test set (Figure.1b). Similar results were also observed in the entire WT Bf-CM set that high risk group had shorter OS than low risk group (47.9 months vs. 170.2 months, p<0.001) and the AUC was 0.780 (Figure.1c). Risk-related profiles revealed that more deaths occurred and expression of signature genes increased with the raise of risk scores (Figure.1d).
3.3 Independence evaluation of the signature in prognostic prediction from clinical and pathological factors
Traditional clinical and pathological factors have been well-clarified to be associated with the prognosis of melanoma patients, including age, gender, pathological stage, Breslow depths and radiation therapy 1,25. Using univariable and multivariable analysis, age and stage were shown to be significantly associated with the prognosis of 211 WT Bf-CM patients in addition to the proposed signature (Figure.2a). Therefore, age and stage were chosen to be compared with the proposed signature for its independence. As is shown in Figure.2b, when applied in each subgroup of entire set divided according to the risk stratification of each factor, the proposed signature revealed considerable performances in prognosis prediction. These results support the independence of this expression-based prognostic signature.
3.4 Comparison of the signature with other known prognostic biomarkers
Considering the burden brought by CM, several studies have proposed other biomarkers for CM prognosis prediction 8,26-30. Although none of them specifically aimed at WT Bf-CM, we compared our signature with them to examine the predictive performances in WT Bf-CM prognosis. As is shown in Figure. 3a and Table. 3, the ROC analysis of time-dependent analysis in WT Bf-CM cohort proved that our signature had significant higher AUC than other biomarkers, and might be superior in WT Bf-CM prognosis prediction exclusively.
3.5 Construction and evaluation of a prognostic nomogram for WT Bf-CM
Nomogram has been a welcomed quantitative predictive tool for clinical practice. Based on the results from univariable and multivariable analysis above (Figure.2a), a nomogram to predict 3−year and 5-year survival of WT Bf-CM patients was constructed based on the risk scores of the proposed signature, age and pathological stage (Figure. 3b). The C-index of this nomogram was calculated to be 0.7505 and calibration plot was also drawn, exhibiting an acceptable accuracy for prediction (Figure. 3c). Additionally, DCA was performed to assess the efficiency of the proposed nomogram (Figure. 3d, e). When the threshold probability was lower than 25%, our nomogram would be more beneficial and when threshold probability higher than 25%, our signature would be more beneficial. These results supported the promising complementary use of our signature and nomogram.
3.6 Functional analysis based on the signature grouping
Cellular and molecular mechanism accounting for the prognosis of melanoma have always been a hot topic of cancer research. Thus, we aimed to find out what biological process might explain the predictive function of this signature. 865 genes were defined as DEGs between high risk and low risk group (Figure. 4a) and interactions of 527 proteins provided a landscape of potential molecular mechanisms (Figure. 4b). Using Mcode, a six-gene module with highest score (score 6) was selected as core module, which was closely associated with mRNA processing, splicing and metabolic process (Figure. 4b). BP enrichment analysis showed various processes were involved, including response to oxygen levels, RNA splicing, regulation of developmental growth, etc. (Figure. 4c). Besides, GSEA analysis indicated that the high-risk group was associated with the up-regulation of GO_NUCLEAR_SPECK (NES = 1.361, p = 0.026), GO_PYRIMIDINE_NUCLEOTIDE_BIOSYNTHETIC_PROCESS (NES = 1.573, p = 0.039), REACTOME_MITOTIC_PROMETAPHASE (NES = 1.413, p = 0.041) and TATA_01 (NES = 1.589, p = 0.005) (Figure. 4d). Taken together, these results supported that the predictive function of this signature might be related to the up-regulation of RNA splicing, transcription, and cellular proliferation in the high-risk group, and upregulation of these processes have been demonstrated to be linked to malignancy of cancer 31,32.
3.7 Evaluation for the therapeutic responses based on signature grouping
The treatment of cutaneous melanoma has made a considerable progress recently, especially for the application of immune-therapy and targeted drugs. Hence, we wondered if the proposed signature could provide clues for the therapeutic responses of WT Bf-CM patients. By estimating the immune infiltration of the WT Bf-CM cohort, we found that little differences in 24 sorts of immune cells infiltration (Figure. 5a). Similarly, the comparison of the expression of 9 immune-markers showed little differences between high- and low- risk groups (Figure. 5b). Of note, although the checkpoint therapies targeting programmed cell death 1 (PD-1), programmed cell death-ligand 1 (PD-L1) and cytotoxic T lymphocyte antigen 4 (CTLA-4) have been recommended for the treatment of WT Bf-CM, their expression didn’t seem to be high both in high- or low-risk groups, suggesting the poor responses of checkpoint therapy in WT Bf-CM patients. Interestingly, the analysis of responses to chemo-therapies showed significant differences of 9 chemotherapeutic drugs in estimated IC50 between high and low-risk groups, along with an increased sensitivity to these 9 chemotherapies in high risk group (Figure. 5c). Among these 9 chemotherapeutic drugs, the targets of AKT.inhibitor.VIII, AS601245, GDC0941, MK.2206, PF.02341066 (Crizotinib), Rapamycin and Temsirolimus are involved in the regulation of PI3K/Akt/mTOR pathway 33-35, and this axis has known to play a critical role in cell proliferation and RNA processing 36,37, further supporting the results from functional analysis. Altogether, these findings support that this signature could predict the responses of PI3K/Akt/mTOR pathway based drugs in WT Bf-CM patients, besides predicting their prognosis.