Identification of Immune-Related LncRNA for Predicting Survival in Skin Cutaneous Melanoma

Qi Jiayi1, Huang XuMiao2, Lin Wei3, Ao YaHui3, Zhao Wenting2* and Cao Chunhao3*

1Chongqing Medical University, Chongqing, 400016, China

2Hubei University of Chinese Medicine

3Department of Integrated Traditional Chinese and Western Medicine, The First Affiliate Hospital of Chongqing Medical University, Chongqing, 400016, China

Qi Jiayi and Huang XuMiao contributed equally to this work.

ABSTRACT

Immune-related lncRNAs (LncRNA) play an essential role in tumorigenesis and progression, and sought to develop an immune-related lncRNA signature to assess the value of immune-related lncRNAs to discover their prognostic values. RNA-seq and clinical information were obtained from the Cancer Genome Atlas (TCGA) database, which aims to screen for immune-associated lncRNAs based on constructing a risk scoring formula. Besides, in virtue of univariate and multivariate Cox regression, immune-associated lncRNA signatures formed its structure. In addition, in virtue of gene ontology (GO) and KEGG pathway enrichment analysis genes were given functional identification. In total, there were 102 differentially expressed immune-related lncRNAs which were identified, and there were four genes (USP30-AS1, U62317.1, ZEB1-AS1 and MIR205HR) which were screened from the gene signature to form a prognostic signature model. Based on the survival analysis results, patients with high-risk scores had poor survival outcomes. To conclude, we constructed a predictive signature of four immune-related lncRNAs, which can provide assistance for prognosis prediction and immunotherapy of SKCM.

Article Information

Received 20 June 2022

Revised 28 June 2022

Accepted 09 July 2022

Available online 30 September 2022

(early access)

Published 24 November 2023

Authors’ Contribution

Skin cutaneous melanoma, Long non-coding RNA, Risk model, The Cancer Genome Atlas (TCGA)

Key words

Coilia nasus, SSR markers, Transcriptome.

DOI: https://dx.doi.org/10.17582/journal.pjz/20220620210604

* Corresponding author: fanquli88@gmail.com, 1932709377@stmail.hbtcm.edu.cn

030-9923/2024/0001-01 $ 9.00/0

Copyright 2024 by the authors. Licensee Zoological Society of Pakistan.

This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Introduction

Skin cutaneous melanoma (SKCM), which originates from melanocytes in the epidermis, is the most common subtype of melanoma (Schadendorf et al., 2018; Leeneman et al., 2021). Most melanomas that are limited to the primary site have a good outcome. However, the median survival for malignant melanoma is only six months, accounting for 75% of skin cancer deaths (Arslanbaeva and Santoro, 2020; Kumar et al., 2021). Because malignant melanomas exhibit very diverse clinicopathological and cytological changes, it can be difficult to identify some unusual malignant melanomas based on pathological variants alone (Cabrera and Recule, 2018). Early identification of aggressive melanoma is essential for developing treatment strategies and predicting disease progression (Eisenstein et al., 2018).

Long non-coding RNA (lncRNA) is a type of RNA without an open reading frame and does not encode proteins, which transcript length is more than 200 nucleotides (Wang et al., 2020). By interacting with proteins, RNA, and DNA in different cellular environments and biological processes, lncRNAs can regulate the expression of nearby genes or genes in other parts of the cell. Recently, research in oncology has focused more on the role of LncRNA in various types of cancer, including SKCM. For instance, CASC2, LINC00961, LINC00459, and MEG3 have been shown as tumor-suppressor lncRNAs in melanoma cancer, while SLNCR1, ANRIL, and SAMMSON exert their oncogenic function (Chen et al., 2020).

Since the incidence of cutaneous melanoma is increasingly growing, it has become more critical to predict the outcome of the disease. Therefore, we used the TCGA dataset to develop a prognostic model for SKCM. Moreover, we explored the relationship between the screened LncRNA and SKCM, which can analyze its diagnostic effectiveness, chemotherapy efficacy, and tumor immune infiltration.

MATERIALS AND METHODS

TCGA dataset and clinical information of patients with SKCM

Gene expression quantification data and corresponding clinical information were acquired from The Cancer Genome Atlas (Li et al., 2021) (https://portal.gdc.cancer.gov/), such as age, gender, tumor stage, survival information, etc. Import the downloaded data into the Perl script for sorting and conversion, and got the gene expression matrix of the above sample.

Identifying immune-related and critical lncRNA

We obtained the immune-related genes from Gene Set Enrichment Analysis (GSEA) website (Yu et al., 2020). The overlapping genes of TCGA data set and immune-related lncRNAs are screened by immune gene-lncRNA co-expression method and correlation analysis in the R software (version 3.6.2).

Prognostic immune-related lncRNAs screening

The differential gene expression data and patient survival data were combined into an expression matrix by R language, and duplicate samples in the obtained expression matrix were excluded. Then a univariate Cox regression analysis was performed by the “survival” package in R language with a screening criterion of P <0.01.

Establishment of immune-related lncRNA signature model

The first 30 lncRNAs selected were analyzed by a multivariate Cox model, and the risk model was constructed with risk value = ∑ (gene coefficient × gene expression amount). Based on the median risk value, patients were divided into high-risk groups and low-risk groups (Sun et al., 2021).

Validation of the risk score model and establishment of the signature for prognosis prediction

Kaplan-Meier curve analysis to estimate the survival difference between high-risk and low-risk groups. R package survival ROC was used to assess the sensitivity and specificity of the lncRNA prognostic model.

Principal components analysis

Principal Components Analysis (PCA) was used to analyze the different Inc-RNA expression patterns between low-risk and high-risk groups (Groth et al., 2013).

Gene set enrichment analysis

GSEA is a calculation method that assesses statistical differences in genes between groups (Zhang and Huang, 2020). The enrichment map is used to visualize GSEA.

GO and KEGG enrichment analyses

To validate the biological functions of the screened lncRNAs. We used R software’s ClusterProfiler package to conduct GO enrichment and KEGG enrichment analyses.

RESULTS

Identification of differentially expressed lncRNAs

Fragments per kilobase per million (FPKM) normalized expression used as the RNA-seq results of 471 SKCM samples and one normal sample obtained from TCGA. One hundred two differentially expressed lncRNAs were identified by median mapping, of which the first 31 are displayed in Supplementary Figure S1. Then, all the above immune-related lncRNAs were analyzed by multivariate Cox regression analysis (P <0.05); 4 immune-related LncRNAs risk evaluation models were established. (Fig. 1).

 

Validation of the 4 immunes-related lncRNAs signature for survival prediction

We calculated the risk score for each sample based on the expression of 4 IRlncRNAs: risk score= (-0.487029792* expression value of USP30-AS1)+(0.465302409 expression value of MIR205HG)+(-0.709815123* expression value of U62317.1)+(-0.378677217* expression value of ZEB1-AS1). Optimal cutoff values for risk scores were calculated using the survminer software package. Patients were divided into high-risk and low-risk groups. Figure 2A, B demonstrate the risk score and survival rate for each case. The heatmap shows that MIR205HG had the highest expression in the high-risk group, while the remaining nine lncRNAs had moderate to high expression levels in the low-risk group (Fig. 2C). Kaplan-Meier plots showed that patients with high-risk scores had a worse probability of OS (Fig. 2D).

 

Clinical evaluation by risk assessment model

We used univariable and multivariable Cox proportion hazard regression models for clinical evaluation of risk assessment models. The analysis showed that age (HR= 1.012, P < 0.05), T-stage (HR= 1.433, P < 0.001), N-stage (HR= 1.728, P < 0.001), and risk score (HR= 1.459, P < 0.001) were significantly associated with OS in the univariable COX regression analysis (Fig. 3A). Similarly, risk group (HR= 1.592, P < 0.001), age (HR= 1.012, P < 0.05), grade (HR= 1.602, P < 0.001) T-stage (HR= 1.433, P < 0.001) and N-stage (HR= 1.728, P < 0.001) maintained their prognostic values in multivariate stepwise cox regression analysis (Fig. 3B). Also, we compared the 5-year ROC curves with other clinical characteristics (Fig. 3C). The results showed that the area under the ROC (receiver operating characteristic) curve (AUC) was greater than 0.5. We further analyzed the expression of the first 6 candidate lncRNAs in relation to different clinicopathological factors (eg. T, N, M, stage 7 AJCC staging, molecular typing, etc.). Our results confirmed that AP003392.1 (p < 0.01), AL022316.1 (P< 0.001), LINC01315 (p < 0.05), and AL161785.1 (p < 0.001) were differentially expressed between different T grades (Fig. 3D).

 

PCA analysis

Information on a total of 471 patients was obtained from the TCGA database, and samples with survival less than 30 days were excluded and included in the PCA analysis. Figure 4A shows the PCA results for all 56753 genes, and Figure 4B shows the PCA results for 10 immune-related genes. The high-risk and low-risk populations are indicated by red and green dots, respectively. PCA analysis of the samples from the TCGA dataset showed significant differences between the samples after RS clustering before and after correction.

 

GSEA for functional annotation of the immune-related risk signature

Gene set enrichment analysis for functional annotation was performed using GSEA software. The results showed a further enhancement of immune response and immune effector processes in the low-risk population compared to the high-risk population (Fig. 5).

GO and KEGG enrichment analyses for the differentially expressed genes between the high-risk group and the low-risk group

In order to further explore the potential functions of the prognostic model, we performed differential expression mRNA analysis on the high-risk and low-risk groups. Biological processes (BP), cell components (CC), and molecular functions (MF) all reveal that they are closely related to tumor immunity.

The top three GO terms for the biological processes were MHC class II protein complex binding, chemokine receptor binding, and cytokine binding. The top three GO terms for the cellular components were clathrin-coated vesicle, coated vesicle membrane, and endolysosome membrane. The top three GO terms for molecular functions were endolysosome membrane, lymphocyte differentiation, and positive regulation of interferon-gamma production (Fig. 6A).

 

KEGG pathway analysis confirmed that primary immunodeficiency, hematopoietic cell lineage and cytokine-cytokine receptor interaction are the most significant enrichment pathways (Fig. 6B). KEGG pathway analysis confirmed that primary immunodeficiency, Hematopoietic cell lineage and cytokine-cytokine receptor interaction are the most significant enrichment pathways (Fig. 6B).

DISCUSSION

At present, it is universally acknowledged that the diagnosed melanoma has been on the increase, which might result from the number of people who have to work outdoors has been growing. That is, being exposed to excessive ultraviolet radiation. With the application of the introduction of targeted therapies and immunotherapy, the survival rate of diagnosed patients has witnessed steady improvement (Pan et al., 2021; Pielesz et al., 2020; Ben Daya et al., 2021). However, treatment outcomes for the patients with the advanced disease remain poor. Over the past decade, the value of molecular biomarkers in predicting disease outcomes, which has become increasingly significant, has attracted the attention of quite a number of investigators. According to the research the available, S100B and LDH are expressing striking serological biomarkers associated with melanoma in predicting disease progression to advanced stages (Eisenstein et al., 2018). With these introduced findings, these trends are emphasizing the importance of identifying biomarkers for prognosis and predicting treatment response. In addition, despite that these patients have similar clinical risk factors molecular tumor characteristics is also able to exert effects upon treatment response and survival in SKCM patients. It eventually leads to the inappropriateness of the traditional clinical-histological classification of T- N- and M-stages, and on that basis, we should further search for prognostically relevant molecular markers.

 

It has been found that LncRNAs can exert a key role in tumorigenesis in various cancers. Additionally, immune-related RNAs are actively involved in regulating the immune status and tumor microenvironment, which may affect prognosis and therapy. Many researchers, who have explored the relationship between immune-related lncRNAs and tumors, have developed risk models on that basis. For example, Liu et al. (2020) who modeled the risk of immune-related lncRNAs in pancreatic cancer, implemented indepth investigation on their part in predicting pancreatic cancer prognosis and immune landscape. It has been exemplified that the lung adenocarcinoma risk model developed by Li et al. (2020) is functioning as valid instrument for risk stratification and individual prognosis assessment. Thereby, it is of a crucial reference for individualized clinical treatment of patients (Xin et al., 2021).

With analyses on the expression profile of cutaneous melanoma patients in the TCGA database, 102 prognostics differentially expressed lncRNAs were screened by co-expression analysis. Univariate Cox regression analysis and multivariate Cox proportional risk analysis was committed to determining the prognostic signatures of four immune-associated lncRNAs (USP30-AS1, U62317.1, ZEB1-AS1, and MIR205HR). Subsequently, K-M curve analysis was performed on patients in the high-risk and low-risk groups. Accordingly, significant differences were found between the high-risk and low-risk groups. In addition, as revealed in the functional enrichment analysis which we performed, these lncRNAs is exerting an indispensable role in the development and progression of SKCM, with multiple biological functional differences between the high- and low-risk groups, including cytokine receptor activity and immune recognition receptor activity. Therefore, we embrace the idea that this model can effectively identify patients at high risk of SKCM.

Among the prognostic signatures established by our group, MIR205HR gives a strikingly high expression in the high-risk group, while USP30-AS1, U62317.1, and ZEb1-AS1 presented high expression in the low-risk group. MIR205HG is a lncRNA containing the miR-205 gene and miR-205 is a miRNA expressed in epithelial cells (Du et al., 2019). Previously, miR-205 conducted differential expression in melanoma tissues from normal tissues. Simultaneously, miR-205 could synergistically affect vascular endothelial growth factors and others with miR-203, miR-9, and miR-15b expression, which aim to regulate self-renewal and epithelial-mesenchymal transition processes. miR-205 levels are up-regulated, which leads to an increase of cell proliferation, migration, and invasion (De Oliveira et al., 2021). It has been demonstrated that MIR205HG levels were more significantly upregulated in melanoma cell lines compared to normal human melanocytes. In virtue of the miR-299-3p/VEGFA axis, the MIR205HG was able to melanoma growth (Guo et al. 2021). High expression of USP30-AS1 specializes striking relevance to poor overall survival prognosis in cervical cancer patients. More than that, via the introduction of miR-299-3p to up-regulate PTP4A1, USP30-AS1 functions as a ceRNA and exacerbates the oncogenic potential of cervical cancer cells both in vitro and in vivo (Chen et al., 2021). However, USP30-AS1 and SKCM does not reveal clear relevance. LncRNA ZEB1-AS1 is up-regulated in melanoma. LncRNA ZEB1-as1 can perform the duty of activating ZEB1 (zinc finger E-box binding homeobox 1) gene expression, thereby affecting the aggressiveness and phenotypic transition of melanoma. It highlights its consistence with our research results. So, it is an indispensable commitment to explore the relationship between LncRNA ZEB1-AS1 and SKCM at a more profound level (Siena et al., 2019). Up to now, no research has been able to unfold the association between U62317.1 and cancer. Despite that, our speculation is that this lncRNA may be involved in the development of SKCM, which accounts for more necessary research for testing this hypothesis. The current manuscript has some limitations as follows. In this study we have evaluated the role of lncRNAs by using bioinformatics analysis and this means that the data represented in this study might not impact on clinical and experimental studies. Nevertheless, we should in the meanwhile also be aware of the shortcomings of this signature; thereby, more clinical data are need to confirm the reliability of this result, and more detailed segmentation is demanded clinical practice, rather than just divide it into high-risk and low-risk groups. On the other hand, due to limitation of resources we could not analyze the subtypes of melanoma.

CONCLUSION

In this study, we established a prognostic signature composed of four immune-associated lncRNAs. This risk model, which can contribute new associations between SKCM and immune-related lncRNAs, can meanwhile provide new strategies for immunotherapy of SKCM.

Supplementary material

There is supplementary material associated with this article. Access the material online at: https://dx.doi.org/10.17582/journal.pjz/20220620210604

Statement of conflict of interest

The authors have declared no conflict of interest.

REFERENCES

Ao, H., Li, H.W., Lu, L.K., Fu, J.X., Han, M.H., Guo, Y.F., and Wang, X.T., 2021. Sensitive tumor cell line for annonaceous acetogenins and high therapeutic efficacy at a low dose for choriocarcinoma therapy. J. Biomed. Nanotechnol., 17: 2062-2070. https://doi.org/10.1166/jbn.2021.3175

Arslanbaeva, L.R., and Santoro, M.M., 2020. Adaptive redox homeostasis in cutaneous melanoma. Redox Biol., 37: 101753. https://doi.org/10.1016/j.redox.2020.101753

Ben-Daya, S.M., Paul, V., Awad, N.S., Al-Sawaftah, N.M., Al-Sayah, M.H., and Husseini, G.A., 2021. Targeting breast cancer using hyaluronic acid-conjugated liposomes triggered with ultrasound. J. Biomed. Nanotechnol., 17: 90-99. https://doi.org/10.1166/jbn.2021.3012

Cabrera, R., and Recule, F., 2018. Unusual clinical presentations of malignant melanoma: a review of clinical and histologic features with special emphasis on dermatoscopic findings. Am. J. clin. Dermatol., 19: 15-23. https://doi.org/10.1007/s40257-018-0373-6

Chen, C.L., Tsai, C.Y., Chen, Y.S., Lin, T.Y., Hsu, Y.J., Wei, S.Y., and Chen, Y.C., 2020. Two-stage patterned cell-based treatments for skin regeneration. J. Biomed. Nanotechnol., 16: 1740-1754. https://doi.org/10.1166/jbn.2020.3003

Chen, M., Chi, Y., Chen, H., and Zhao, L., 2021. Long non-coding RNA USP30-AS1 aggravates the malignant progression of cervical cancer by sequestering microRNA-299-3p and thereby overexpressing PTP4A1. Oncol. Lett., 22: 1-12. https://doi.org/10.3892/ol.2021.12766

De Oliveira, T.D., Travassos, L.R., Arruda, D.C., and Tada, D.B., 2021. Intracellular targeting of poly lactic-co-glycolic acid nanoparticles by surface functionalization with peptides. J. Biomed. Nanotechnol., 17: 1320-1329. https://doi.org/10.1166/jbn.2021.3108

Du, Q., Hoover, A.R., Dozmorov, I., Raj, P., Khan, S., Molina, E., and van Oers, N.S., 2019. MIR205HG is a long noncoding RNA that regulates growth hormone and prolactin production in the anterior pituitary. Dev. Cell., 49: 618-631. https://doi.org/10.1016/j.devcel.2019.03.012

Eisenstein, A., Gonzalez, E. C., Raghunathan, R., Xu, X., Wu, M., McLean, E.O., and Alani, R.M., 2018. Emerging biomarkers in cutaneous melanoma. Mol. Diagn. Ther., 22: 203-218. https://doi.org/10.1007/s40291-018-0318-z

Groth, D., Hartmann, S., Klie, S., and Selbig, J., 2013. Principal components analysis. Comput. Toxicol. Methods Mol. Biol., 930: 527. https://doi.org/10.1007/978-1-62703-059-5_22

Hou, Q., Huang, J., Xiong, X., Guo, Y. and Zhang, B., 2022. Role of nutrient-sensing receptor GPRC6A in regulating colonic group 3 innate lymphoid cells and inflamed mucosal healing. J Crohns Colitis, 16: 1293-1305. https://doi.org/10.1093/ecco-jcc/jjac020

Lei, J., Dong, Y., Hou, Q., He, Y., Lai, Y., Liao, C. and Zhang, B., 2022. Intestinal microbiota regulate certain meat quality parameters in chicken. Front. Nutri. (Lausanne), 9: 747705. https://doi.org/10.3389/fnut.2022.747705

Li, H., Peng, R. and Wang, Z., 2018. On a diffusive susceptible-infected-susceptible epidemic model with mass action mechanism and birth-death effect: Analysis, simulations, and comparison with other mechanisms. SIAM J. Appl. Math., 78: 2129-2153. https://doi.org/10.1137/18M1167863

Li, M., Hu, M., Zeng, H., Yang, B., Zhang, Y., Li, Z., and Ming, Y., 2021. Multifunctional zinc oxide/silver bimetallic nanomaterial-loaded nanofibers for enhanced tissue regeneration and wound healing. J. Biomed. Nanotechnol., 17: 1840-1849. https://doi.org/10.1166/jbn.2021.3152

Liu, Z., Su, W., Ao, J., Wang, M., Jiang, Q., He, J. and Ji, M., 2022. Instant diagnosis of gastroscopic biopsy via deep-learned single-shot femtosecond stimulated Raman histology. Nat. Commun., 13: 4050. https://doi.org/10.1038/s41467-022-31339-8

Pan, Y., Hao, Y., Xiao, Y., Shi, K., Qu, Y., and Qian, Z., 2021. Injectable soft tissue nano/micro fillers for facial reconstruction. J. Biomed. Nanotechnol., 17: 1-17. https://doi.org/10.1166/jbn.2021.3011

Pielesz, A., Biniás, D., Machnicka, A., Bobínski, R., Laane, H.M., Tjon, J., and Waksmánska, W., 2020. Stabilized ortho-silicic acid as a modifier of tissue preliminary comparative biomaterial studies. J. Biomed. Nanotechnol., 16: 538-547. https://doi.org/10.1166/jbn.2020.2912

Schadendorf, D., van Akkooi, A.C., Berking, C., Griewank, K.G., Gutzmer, R., Hauschild, A., and Ugurel, S., 2018. Melanoma. Lancet, 392: 971-984. https://doi.org/10.1016/S0140-6736(18)31559-9

Siena, Á.D.D., Plaça, J.R., Araújo, L.F., de Barros, I.I., Peronni, K., Molfetta, G., and Silva, W.A., 2019. Whole transcriptome analysis reveals correlation of long noncoding RNA ZEB1-AS1 with invasive profile in melanoma. Sci. Rep., 9: 1-11. https://doi.org/10.1038/s41598-019-47363-6

Sun, L., Xiang, H., Ge, C., Chen, X., Zhang, Q., Zhang, Y., and Miao, X., 2021. A nanocrystals-based topical drug delivery system with improved dermal penetration and enhanced treatment of skin diseases. J. Biomed. Nanotechnol., 17: 2319-2337. https://doi.org/10.1166/jbn.2021.3202

Wang, Y., Sun, B., Wen, X., Hao, D., Du, D., He, G., and Jiang, X., 2020. The roles of lncRNA in cutaneous squamous cell carcinoma. Front Oncol., 10: 158. https://doi.org/10.3389/fonc.2020.00158

Xin, S., Wu, Y., Huang, Z., Huang, Y., Jia, B., and Zhao, J., 2021. Engineering cell membrane-based nanovesicles for melanoma tumor treatment. J. Biomed. Nanotechnol., 17: 838-845. https://doi.org/10.1166/jbn.2021.3072

Yu, Q., Huang, Y., Zhu, C., Wu, X., Tai, Z., Xie, X., and Zhu, Q., 2020. Combination of microneedles and mf59 adjuvant as a simple approach to enhance transcutaneous immunization. J. Biomed. Nanotechnol., 16: 1776-1786. https://doi.org/10.1166/jbn.2020.3007

Zhang, D., and Huang, Q., 2020. Encapsulation of astragaloside with matrix metalloproteinase-2-responsive hyaluronic acid end-conjugated polyamidoamine dendrimers improves wound healing in diabetes. J. Biomed. Nanotechnol., 16: 1229-1240. https://doi.org/10.1166/jbn.2020.2971