Identification of ERS-related lncRNAs
By constructing the co-expression network of mRNA and lncRNA, a total of 14,142 ERS-related lncRNAs were obtained. Then, 39 ERS-related lncRNAs were found to be associated with survival prognosis in skcm patients using univariate Cox regression analysis (P < 0.001) (Figure.1A). We analyzed the differential expression of 39 ERS-related lncRNAs in normal tissues and SKCM (Figure. 1B), and analyzed the correlation between these lncRNAs and PD-L1 (Figure. 1C).
Consensus clustering of SKCM patients based on ERS-related prognostic lncRNAs
Based on the differential expression of ERS-related prognostic lncRNAs in different SKCM patients, patients were divided into different subgroups using consensus clustering. Based on the similarities shown by ERS-related prognostic lncRNA expression levels, we found that each cluster had the highest stability when k = 2 (Figure. 2A). SKCM patients were divided into cluster 1 of 168 patients and cluster 2 of 290 patients (Fig. 3(a)). Compared with skcm patients in cluster 2, patients in cluster 1 had a higher overall survival rate (Figure. 2B, p≤0.001). Interestingly, we also found some correlation between patients in cluster 2 and advanced clinical stage (Figure. 2C).
Consensus Clustering related to tumor microenvironment
By immune scoring patients in cluster 1 and cluster 2, we found differences between the two clusters in the level of immune cell infiltration, which gave us a better understanding of the role of ERS-related prognostic lncRNAs in the immune microenvironment of SKCM. The scores in cluster 1 were higher than those in cluster 2, whether it wass immune score, matrix score, or estimated score (Figure 3A-C). Figure 3D presents the infiltration abundance of 22 types of immune cells in cluster 1 and cluster 2. We found significant differences in the degree of infiltration of 12 immune cells in cluster 1 and cluster 2, which were B cells memory, Plasma cells, T cells CD8, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, NK cells resting, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting and Mast cells activated(supplement figure1).
Enrichment analysis of 2 SKCM subtypes
We explored the underlying regulatory mechanisms responsible for the differences between the two clusters of SKCM patients by performing KEGG enrichment analysis in GSEA. The top five pathway enrichments significantly associated with cluster 1 and cluster 2 are showed in Figure 5. Chemokines, cytokine receptor interactions, and JAK-STAT signaling pathways were found to be associated with cluster 1, While cluster 2 is mainly associated with RNA polymerase and aminoacyl tRNA biosynthesis.
Construction of a ERS-related prognostic lncRNA signature
The 39 prognostic ERS-related lncRNAs obtained were further analyzed by lasso regression, and a prognostic diagnostic risk model consisting of 10 ERS-related lncRNAs(AATBC, ZEB1-AS1, LINC01871, AC093726.1, AC021188.1, AC092747.4, AC009495.2, AL157871.2, AC242842.1 and AC067930.4) was obtained. The coefficient and partial likelihood deviance of prognostic signature are presented in Figures 6A and Figures 6B. By random assignment, 230 SKCM patients were assigned to the training cohort and 228 SKCM patients were assigned to the testing cohort. The 10 ERS-related lncRNAs and their corresponding risk coefficient are showed in table2. Risk score for each patient: Risk score = [(0.051* Exp (AATBC) + (-0.276* Exp (ZEB1-AS1) + (-0.118* Exp (LINC01871) + (-0.004* Exp (AC093726.1) + (-0.101* Exp (AC021188.1) + (-0.017* Exp(AC092747.4) + (0.142* Exp(AC009495.2) + (-0.045* Exp(AL157871.2) + (-0.067* Exp(AC242842.1)+ (-0.127* Exp(AC067930.4)]. Then, according to the risk score, SKCM patients were divided into high-risk and low risk groups. Figure 7 shows the distribution of survival probability(Figure 7A and Figure 7E), survival time(Figure 7C and Figure 7G), risk score((Figure 7B and Figure 7F) and lncRNA signature ((Figure 7D and Figure 7H)) of each SKCM patient in the training and testing cohort. In both the training cohort and the test cohort, the prognosis was better in the low-risk group than in the high-risk group, both in terms of survival probability and survival time(Figure7A).
Univariate regression analysis showed that T stage, N stage, clinical stage and risk score were associated with survival prognosis in both training cohort and testing cohort(Figure 8B and Figure 8E). Further analysis by multivariate regression showed that T stage, N stage and risk score were still significantly associated with prognosis(Figure 8C and Figure 8F). Therefore, we believed that the risk score can serve as an independent prognostic risk factor for patients with SKCM. The ROC curve analysis showed that the risk score had a high predictive accuracy for patient prognosis in both training cohort and testing cohort, with areas under curve(AUC) of 0.709(Figure 8A) and 0.719(Figure 8B), respectively.
In addition, we validated the prognosis of risk scores in different groups of SKCM patients and found that in both male (Supplementary Figure. 2C) and female (Supplementary Figure. 2D), and SKCM patients older than 65 years or younger (Supplementary Figure. 2A-B) , the high-risk group had a worse prognosis.
Furthermore, among SKCM patients with advanced clinical stage (Supplementary Figure. 3A), pT3-4 (Supplementary Figure. 3B), pM0 (Supplementary Figure. 3C), and pN1-3 stage (Supplementary Figure. 3D), the low-risk group had a better prognosis compared with the high-risk group,.
Correlations between risk scores and clinical characteristics
AATBC and AC009495.2 were highly expressed in the high-risk group, and the remaining 8 lncRNA signature(ZEB1-AS1, LINC01871, AC093726.1, AC021188.1, AC092747.4, AL157871.2, AC242842.1 and AC067930.4) were mainly highly expressed in the low-risk group(Figure 8A). The heatmap in Figure 8 also showed the differential expression of patients in 2 different risk groups in terms of pN stage, clinical stage, immune score, and cluster subtype. The risk score increased when pT stage increased. In addition, groups with higher immuneScore had higher risk scores compared to group with lower immuneScore.
Risk Score Associated with Immune Infiltration
Risk scores were negatively correlated with the infiltrating abundance of B cells memory, Macrophages M1, T cells CD4 memory activated, T cells CD8 and T cells follicular helper(Figure 9A-E). But the risk score was positively correlated with the infiltrating abundance of immune cells such as Macrophages M0, Macrophages M2 and Mast cells resting(Figure 9F-H).
Response of high-risk and low-risk patients to chemotherapy, targeted therapy, and Immunotherapy
The pRRophetic algorithm was used to predict the IC50 of three common chemotherapeutic agents (docetaxel,cisplatin and paclitaxel) and one common targeted therapeutic agents in high-risk and low-risk patients. There were significant differences in the sensitivity to cisplatin, paclitaxel, and sorafenib between patients in the high-risk group and the low-risk group (Figure 9I-K). Figure 9I-K reveals that patients in the low-risk group are more sensitive to the chemotherapy drugs paclitaxel and cisplatin than those in the high-risk group. However, patients in the high-risk group were more sensitive to sorafenib than those in the low-risk group (Figure 9K). To investigate the potential susceptibility of SKCM patients to immune checkpoint (PD-L1, CTLA-4, HAVCR2, LAG-3, PDCD1, or VSIR) inhibitors, we compared the differences in immune checkpoint gene expression between high-risk and low-risk patients. The results showed that the expression levels of these six immune checkpoint genes in the low-risk group were higher than those in the high-risk group (Figure 9 L-Q). These data suggest that low-risk patients may respond better to immune checkpoint inhibitors targeting PD-L1, CTLA-4, HAVCR2, LAG-3, PDCD1 or VSIR than patients in the high-risk group.
Different ERS statuses in low-risk and high-risk groups
By performing principal component analysis (PCA) on risk models of 10 ERS-related lncRNAs, 1613 ERS-related coding genes and genome-wide expression profiles from TCGA (Supplementary Fig. 4), we found that low-risk and high-risk groups were distributed in two different direction (Supplementary Figure 4). It showed that the risk model could divide SKCM patients into two parts, and the ERS status of SKCM patients differed in high-risk and low-risk groups.