Tumor-infiltrating macrophage associated lncRNA signature in cutaneous melanoma: implications for diagnosis, prognosis, and immunotherapy

Macrophage associated lncRNA in melanoma

In TIMER2.0 website (http://timer.comp-genomics.org/), we astonished observed that tumor infiltrated macrophage cells were regarded as a prognostic factor and intimately associated with BRAF mutant in cutaneous melanoma patients (Supplementary Figure 1). Thus, the 3346 macrophage associated lncRNAs in cutaneous melanoma (Supplementary Table 1) were retrieved from the ImmLnc database (http://bio-bigdata.hrbmu.edu.cn/ImmLnc), which affords tools to explore the immune associated function of lncRNAs such as the correlation between lncRNA and immune cell types, lncRNA and pathways and cancer related lncRNAs across 33 types of cancer [19].

Melanomas collection and normal controls

The TCGA-SKCM dataset (https://portal.gdc.cancer.gov) was used to obtain the lncRNAs expression and clinical variables of cutaneous melanoma, and the normal healthy skin tissue samples were obtained from the Genotype-Tissue Expression website (GTEx) (https://www.gtexportal.org/) to match the tumor samples. TPMs values were utilized to describe the expression level of lncRNAs for TCGA-SKCM & GTEx. Five independent melanoma patient datasets were acquired from GEO database for outside verification, which including GSE65904, GSE15605, GSE78220 (Hugo et al. study) [20] and GSE91061 (Riaz et al. study) [21]. Furthermore, to predict immunotherapy response, A prior cohort managed with CTLA-4 and PD-1 inhibitors was obtained from a previously published study [22]. The transcriptome profiles (FPKM) of these datasets were processed by applying their corresponding platform and was transformed into TPMs. Firstly, probes were annotated based on the annotation profile supplied by the platform, and non-matched probes were removed. If several probes were mapped to the same gene, the average value would be used to reflect gene expression. Variations in lncRNAs with low expression levels (TPMs <0.01) will then be eliminated due to their low expression level. Finally, the raw data were log2(× + 1) transformed and quantile normalized. Figure 1 depicts a simplified procedure for the current investigation.

Differential expression analysis of macrophage associated lncRNAs (DEMlncRNAs)

To explore the differentially expressed lncRNAs, samples in TCGA-SKCM & GTEx dataset were divided into melanoma and normal skin groups. The differentially expressed analysis of lncRNAs between melanoma and normal skin was conducted by using ‘limma package’ in R software [23]. The cutoff criterion was the absolute value of log 2-fold change (FC) ≥2 and p-values < 0.05. Besides, the overlapping lncRNAs in differentially expressed lncRNAs and macrophage associated lncRNAs were defined as DEMlncRNAs.

Development of DEMlncRNAs signature using machine learning

To identify the survival-related DEMlncRNAs in melanomas, we randomly categorized the melanoma patients in TCGA-SKCM dataset into train or validation samples with a ratio of 7:3. The survival associated DEMlncRNAs were then selected by using univariate Cox regression methodology in train dataset (p-values < 0.05). Next, multiple machine learnings contained Least Absolute Shrinkage and Selection Operator (LASSO), Random Forests Variable Selection (RF-VS), and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) were applied to identify the essential DEMlncRNAs in the training dataset. LASSO screened the important features based on the misclassification error, SVM-RFE identified key factors by 5-fold cross-validation errors and RF-VS picked out significant variables by out of bag (OOB) error [24]. The common important DEMlncRNAs were selected out for further analysis. Firstly, the diagnostic model was built by neural network algorithm with these identified DEMlncRNAs. To assess the model’s efficiency, the confusion matrix and Receiver operating characteristic curves (ROC) were utilized. The prognostic risk model was then built using multivariate Cox regression analysis. The following is how the risk score is calculated:

$$\text{Risk}\text{score}\text{=}{\displaystyle {\sum }_{i=1}^{\text{N}}(coe{f}_i\times exp{r}_i),}$$

where expri means the DEMlncRNA’s relative expression in the risk model, coefi is stand for the DEMlncRNA’s Cox coefficient, and N represents the number of DEMlncRNA. In the train set, the median value was used as a threshold for the DEMlncRNA signature risk score, and samples were divided into high- and low-risk categories. Kaplan-Meier survival curves and log-rank tests were used to compare the prognoses of the high- and low-risk groups. ROC curves and area under the curve (AUC) measurements were also used to evaluate the model’s 3- and 5-year overall survival prediction accuracy. Furthermore, these DEMlncRNAs signatures were verified in the testing dataset and GSE65904 to establish the robustness of the conclusion.

Pathway enrichment analysis

To investigate the possible involvement of protein-coding genes co-expressed with DEMlncRNAs, the ENCORI website (http://starbase.sysu.edu.cn/index.php) was used. Second, Metascape (http://metascape.org) was utilized to perform functional enrichment of BP and the KEGG terms. In addition, immune-related pathways of DEMlncRNAs were examined in the ImmLnc database (http://bio-bigdata.hrbmu.edu.cn/ImmLnc). Pathways with p-value and adjust p-value < 0.05 were screened out and regarded as significance.

Associations between risk score and other prognostic markers

To calculate the relationships between risk score and other prognostic biomarkers, the correlation of risk score with traditional markers (age, clark level, stage.) and previously established prognostic markers included transcriptomic classification (named “immune”, “keratin” and “MITF-low”), mutation subtype (named “BRAF subtype”, “RAS subtype”, “NF1 subtype” and “Triple Wild-Type subtype”) and BRAF mutant were performed [25, 26]. Moreover, uni- and multi-variate Cox regression were performed in three datasets to analyze the prognostic significance of the risk score and traditional clinical factors to evaluate if the risk score is an independent prognostic factor. Next, to explore the association between risk score and immunological microenvironment, the Immuno-Oncology Biological Research (IOBR) approach was used. The CIBERSORT method was initially used to determine the relative proportions of 22 different kinds of infiltrated immune cells in each patient. Only immune cells have 50% value and patients with P < 0.05 were considered eligible for further analysis. Afterwards, previous user-built signatures associated with immune microenvironment were evaluated by ssGSEA, Principal Component Analysis (PCA), and Z-score methods. Finally, the associations for risk score with immune cell infiltration and immune-associated markers were studied further.

Statistical analysis

Every statistical analyses were carried out by R (version 3.6.0) or Python (version 3.8.0) and installed packages. The “glmnet” package performed the LASSO analysis. The “e1017” package was used to carry out the SVM-RFE approach. The “varSelRF” package was applied to the RF-VS algorithm. The “CIBERSORT” program estimated the CIBERSORT deconvolution. The “survival” and “survivalROC” packages were used to perform Kaplan-Meier (KM) survival analyses. The “IOBR” package estimated the IOBR algorithm. Python’s “sklearn” module was used to run the neural network method. The Spearman algorithm was used to evaluate the correlation test. When comparing two groups, a Wilcoxon test was employed; when comparing more than two groups, Kruskal-Wallis tests was utilized. To investigate the relationship between subgroup and clinicopathological variables, a Chi-square test was performed. The Cox regression analysis yielded hazard ratios (HR) and 95 percent confidence intervals (CI). In all tests, P < 0.05 was used to denote statistical significance.

Availability of data and materials

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Differentially expressed lncRNA in melanoma

The normalized matrix data of lncRNA for cutaneous melanoma and normal samples (not sun exposed) were obtained from the TCGA-SKCM & GTEx, which included 466 melanomas and 233 normal samples. According to the selection criterion, 1640 differentially expressed lncRNAs were screened out in TCGA-SKCM & GTEx dataset, where 650 lncRNAs were significant up-regulation and 990 lncRNAs were significant down-regulation (Supplementary Figure 2A).

Identification of validation of prognostic DEMlncRNAs signature

According the results of overlap in macrophage associated lncRNAs and differentially expressed lncRNAs (Supplementary Figure 2B), 447 DEMlncRNAs were selected for further analysis (Supplementary Table 2). The heatmap of 447 DEMlncRNAs was shown in Supplementary Figure 2C. After excluding the patients without survival information, the TCGA-SKCM dataset were randomly separated into two independent datasets: train samples (n = 313) and validation samples (n = 134). Comparison of the clinical information between train and validation samples, no statistically significant differences were observed in two datasets (Table 1). Next, the relationships between 447 DEMlncRNAs and overall survival time in train sample were assessed by univariate Cox regression analysis, and we discovered 115 DEMlncRNAs were significantly correlated with survival (Figure 2A). In order to find the optimum signatures, machine learning methods were firstly performed to screen out the most important DEMlncRNAs in the train dataset. Combined the feature selection results from LASSO algorithm (Figure 2B), RF-VS algorithm (Figure 2C) and SVM-RFE algorithm (Figure 2D), we observed that six overlapping DEMlncRNAs were shared in these machine learning methods (Figure 2E and Supplementary Table 3). Further comparison analysis between 58 melanomas and 16 normal skin samples in the GSE15605 dataset showed that five overlapping DEMlncRNAs had significant differences (Figure 2F). Furthermore, we used neural network algorithm to well establish a diagnostic model with six DEMlncRNAs in TCGA-SKCM & GTEx dataset (randomly split into train and test dataset with ratio = 7:3), which suggested good accuracies and AUC values in train dataset (accuracy = 0.99 and AUC = 0.99) (Figure 3A, 3D) and test dataset (accuracy = 0.99 and AUC = 1.00) (Figure 3B, 3E). Even in GSE15605 dataset, our diagnostic model still revealed an excellent performance to distinguish normal and tumor patients with accuracy = 0.85 and AUC = 0.93 (Figure 3C, 3F).

Afterwards, the six DEMlncRNAs were subsequently performed to develop a risk model in train dataset by applying multivariate Cox analysis. Each patient’s risk score was finally computed using the risk score formula. The risk model consisted of risk scores distribution, vital status, overall survival (OS) time, and heat map of DEMlncRNAs expression, which respectively illustrated in train (Supplementary Figure 3A), validation (Supplementary Figure 3B) and GSE65904 (Supplementary Figure 3C) datasets. Next, we divided the melanoma patients in the train dataset into high-risk and low-risk groups based on the median value of the risk score. With a significant log-rank test of p < 0.001, the KM survival curves revealed that patients of high-risk group had substantially worse survival than low-risk group (Figure 4A). The ROC curve was drawn to estimate the ability of prediction power, and the AUCs of the curves were 0.619 at 3-year and 0.676 at 5-year respectively (Figure 4B). Moreover, to confirm the robustness of the six DEMlncRNAs based model, it was further verified in validation and GSE65904 datasets via the same cutoff value. The patients in validation and GSE65904 datasets were accordingly divided into high- and low-risk subgroups. KM curves indicated that significantly distinct survival outcome was also observed at high- and low-risk group in validation dataset with log-rank p < 0.001 (Figure 4C) and GSE65904 dataset with log-rank p = 0.026 (Figure 4E). The 3-, 5-year of AUCs were 0.778, 0.764 in validation dataset (Figure 4D) and 0.636, 0.622 (Figure 4F) in GSE65904 dataset respectively.

Independence of the risk score from traditional clinical variables

In uni- and multi-variate Cox regression analysis, the risk score and conventional clinical parameters were utilized to establish if the DEMlncRNAs signature is an independent predictive predictor of OS in various datasets. The train and validation datasets demonstrated that the DEMlncRNAs signature risk score and tumor stage were substantially linked with OS in both uni- and multi-variate Cox analysis (Figure 5A, 5B). In same analyses of an external GSE65904 dataset, the risk score was remained strongly linked with OS (Figure 5C). These findings suggest that the risk score of the DEMlncRNAs signature is a significant predictive factor that is independent of other conventional clinical factors.

Pathway enrichment analysis

Overall, 53 lncRNA-paired protein-coding genes were predicted on the ENCORI website and utilized in BP term and KEGG pathway enrichment analyses. The KEGG enrichment results showed that these paired protein-coding genes were considerably enriched in biological processes linked to RNA transport, the mRNA surveillance pathway, and ribosome synthesis in eukaryotes, transcriptional misregulation in cancer, Herpes simplex infection and Spliceosome (Figure 4G). The BP findings indicated that these paired protein-coding genes were considerably enriched in pathways such as viral gene expression, gene expression regulation, epigenetics, positive control of translation, and 3′-UTR-mediated mRNA destabilization, regulation of telomere maintenance (Figure 4H). Furthermore, Based on ImmLnc database, these DEMlncRNAs significantly associated with immune pathways, including TNF family members, natural killer cell cytotoxicity, interleukins receptor and cytokine Receptors (Figure 4I). The detail results of pathway enrichment analysis were listed in Supplementary Table 4.

Associations between risk score and conventional variables as well as the immunological microenvironment

Box plots demonstrated that clark level, tumor size, radiation, vital status, and metastatic status were associated with risk score (Figure 6A). To investigate the relationships between risk score and immunological microenvironment, the CIBERSORT algorithm was employed to calculate the relative proportions of 22 immune cells. Only 133 melanoma samples and 11 immune cells were selected out for further research based on the selection criteria. In overall, the 11 immune cell subsets in melanoma included memory B cells, naïve B cells, resting Mast cells, and M0, M1, M2 macrophages, resting NK cells, CD 8 T cells, follicular helper T cells, regulatory Tregs T cells and CD4 memory T cells have been activated, whose sum of mean proportions was more than 75% in all melanoma samples (Figure 6B). The differences in immune cell infiltrations between high- and low-risk groups were explored, and the findings revealed that Macrophages M0, M1, CD8, follicular helper, regulatory Tregs and CD4 memory activated T cells were differently infiltrated (Figure 6C). Afterwards, a total of 25 promising immune associated signatures were scored and classified into six categories which contained immune related biomarker, immune microenvironment, immune suppression, immune exclusion, immune exhaustion, T-cell or B-cell receptor (TCR_BCR). The heatmaps delineated a different immune expressed pattern between high and low risk phenotype regardless by using PCA (Figure 7A), ssGSEA (Figure 7B) or Z-score (Figure 7C) methods.

Comparison of risk score and previously identified prognostic indicators

In order to correlate our assessment of risk score with other established prognostic markers, several well established prognostic markers such as transcriptomic classification, mutation subtype and BRAF mutant for cutaneous melanoma were obtained from the previous TCGA study [25]. The TCGA research proven that patient in “immune” subtype of transcriptomic classifications (Figure 8A), “BRAF Hotspot Mutants” of mutation subtypes (Figure 8C) and BRAF mutant group (Figure 8E) have a good prognosis, respectively. The risk distribution of these well established prognostic markers suggested that “MITF low” subtype (Figure 8B), “Triple WT” mutant subtype (Figure 8D) and BRAF wildtype group (Figure 8F) have a higher risk score than “immune” subtype, “BRAF Hotspot Mutants” subtypes and BRAF mutant group, respectively. Most importantly, compared with the 5 years AUC values of these established prognostic markers (transcriptomic classification, mutation subtype, tumor stage and BRAF mutant), our signature can achieve higher accuracy value (Figure 8G).

Potential indicator for melanoma immunotherapy

The relationship between risk score and immune checkpoint gene expression (PD-1 and CTLA-4) was investigated further to assess the potential responsiveness to immunotherapy. According to the Spearman tests, the estimated risk score was substantially negatively linked with PD-1 expression (r = −0.4883; p < 0.0001) (Figure 9A) and CTLA-4 expression (r = −0.2574; p < 0.0001) (Figure 9C). The stratified analysis indicated that PD-1 (Figure 9B) and CTLA-4 (Figure 9D) in the low risk subgroup had significant higher expression levels than those in high risk subgroup. The cross-talk influences of the risk score and immune checkpoint genes on patient survival were further analyzed. Based on the combination of risk score and immune checkpoint genes, these patients were classified into four categories, and KM curve analyses were performed to assess the various survival outcomes among the four subgroups. When compared to the data, the risk score was able to discriminate patients’ outcomes with opposing levels of immune checkpoint genes (Figure 9E, 9F). Patients with high risk scores and low levels of immune checkpoint genes fared the poorest. Patients with low risk scores and high levels of immune checkpoint genes, on the other hand, are more likely to survive longer among the four groups.

Because of the strong association between the risk score and immune checkpoint genes, we postulated that the risk score may be used to predict response to melanoma immunotherapy. As a consequence, we conducted the TIDE method and subclass mapping to predict melanoma who responded to immune-checkpoint inhibitors (CTLA-4 and PD-1) [22]. Strikingly, the low risk group reacts better to anti-PD-1 therapy (Bonferroni adjusted P = 0.007) (Figure 10A). Patients in the high-risk category, on the other hand, do not respond to anti-CTLA-4 therapy (Bonferroni corrected P = 0.039) (Figure 10A). To validate our hypothesis, participants in the Riaz et al. and Hugo et al. melanoma studies who received anti-PD-1 treatment were divided into low and high risk score groups [20, 21]. Notably, individuals in high risk group have a poor prognosis regardless in Hugo et al. cohort (Figure 10B) or Riaz et al. cohort (Figure 10C). Furthermore, in the Hugo et al cohort (Figure 10D) and Riaz et al. cohort (Figure 10E), the non-respond group had a higher risk score than the response group for anti-PD-1 immunotherapy.

Gene set enrichment analysis (GSEA)

GSEA was used to examine the distinct activation of cancer hallmarks between low and high risk groups for investigating the potential biological and molecular mechanism of the lncRNA signature. We found that 10 positive linked pathways were enriched in the low risk group, which included interferon alpha and gamma response, allograft rejection, inflammatory response, epithelial mesenchymal transition, TNFA signaling through NFKB (Figure 10F). Oxidative phosphorylation, MYC targets, estrogen response late, estrogen response early, and KRAS signaling were all actively related with the high risk group (Figure 10F).