Genome-wide screening and immune landscape suggest a potential-m6A-related lncRNA risk signature for predicting prognosis of melanoma

Introduction

Melanoma is the most dangerous form of skin cancer because of its high metastatic potential (1). According to global cancer statistics, morbidity and mortality rates associated with melanoma are increasing yearly (2,3). Based on clinical stage, primary treatment includes surgery, chemotherapy, and molecular targeted therapy. However, for patients with advanced-stage diseases, the therapeutic effect remains poor because of resistance to chemotherapy, aggressive clinical behavior, and high metastatic potential (4). Therefore, novel biomarkers are urgently needed in melanoma for predicting disease progression and prognosis, and to guide therapeutic decisions.

Many researchers have recently confirmed that the immune status plays a vital role in tumorigenesis and cancer progression (5-7). Thus, immunotherapy has been receiving increased attention in cancer treatment. For melanoma, immune checkpoint inhibitors (ICIs) targeting programmed death-ligand 1 (PD-L1) are widely used in clinics (8), especially for patients with advanced melanoma.

Long noncoding RNAs (lncRNAs) do not encode proteins; however, they influence gene expression at different levels. Researchers have confirmed that lncRNAs play a critical role in patients with melanoma. For instance, Wang et al. found that the lncRNA TTN-AS1 promotes melanoma oncogenesis and metastasis by regulating the expression of TTN (9). Xie et al. proved that the lncRNA CAR10 accelerates growth and metastasis by regulating the miR-125b-5p/RAB3D axis in melanoma (10).

N6-methyladenosine (m6A) mediates a major part of RNA methylation and plays a crucial role in RNA splicing, export, and translation, further affecting the occurrence and development of tumors (11,12). m6A modification is controlled by m6A regulators, including “writers,” “readers,” and “erasers” (13). Many studies have shown that m6A methylation affects the development and prognosis of melanoma. Dahal et al. found that METTL3 is upregulated in melanoma tissues and influences tumor invasion/migration through MMP2 (14). Yang et al. found that FTO facilitates tumorigenesis and anti–PD-1 resistance in melanoma (15). However, only a few scientists have explored the potential mechanisms whereby m6A modifications regulate lncRNA-dependent melanoma tumorigenesis and progression.

We thoroughly explored the relationship between potential-m6A-related lncRNAs (pMRlncRNAs), PD-L1 expression level, and tumor immune infiltration in melanoma based on gene expression data and related clinical characteristics data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases in the current study. Furthermore, clustering subgroups and a 10-pMRlncRNA risk signature were developed for prognostic risk stratification and to facilitate decision-making for therapeutic approaches in patients with melanoma. The validity and reliability of the clustering subgroups and pMRlncRNA risk signature were examined using comprehensive bioinformatics analysis. A nomogram plot was constructed for better clinical application.

We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-4402/rc).

Methods

Public dataset and m6A-related genes

The gene expression profiles and related clinical characteristics data of patients with melanoma were downloaded from TCGA database. We used the GSE65904 dataset from the GEO database as the validation cohort. Only melanoma patients with a follow-up time of more than 1 month were included in our study. Patients with incomplete clinical characteristics data were also excluded. Based on some important studies in the field of m6a (16-21), 23 genes act as readers, writers, and erasers are considered to be m6a-related genes (Table S1), and we generated the expression matrixes of these genes. R 3.6.2 (R Core Team) and Perl (Perl Foundation) software were used for data preprocessing. The “limma” (22) and “sva” (23) packages were used to perform data normalization by averaging multiple expressions of the same gene and log2 transformation. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Annotation of pMRlncRNAs

We annotated lncRNAs based on the human reference genome website (GRCh38.p12; https://www.ncbi.nlm.nih.gov/genome). Then, Pearson correlation analysis was implemented to identify pMRlncRNAs with the criteria of P<0.001 and |R| >0.4. Furthermore, univariate Cox regression analysis was implemented to identify pMRlncRNAs related to overall survival (OS) based on TCGA survival data with P<0.05. pMRlncRNAs were regarded as deleterious or protective according to their hazard ratio (HR) values.

Subtype analysis

To explore the biological function of pMRlncRNAs in melanoma, the R package “ConsensusClusterPlus” (1,000 iterations and resample rate of 80%) was used to divide the patients with melanoma into different subtypes in TCGA cohort. A total of 447 patients with melanoma were divided into 2 clusters. Then, the survival and clinicopathological factors of the 2 clusters were analyzed using the R packages “survival” and “pheatmap”. To explore the relationship between PD-L1 and pMRlncRNAs, we compared the expression levels of PD-L1 in the 2 clusters.

Analysis of immune infiltration in the 2 clusters

With the help of the R package “estimate”, we evaluated immune infiltration in patients with melanoma from TCGA cohort. The proportion of 22 immune cell types in melanoma was obtained by CIBERSORT (https://cibersort.stanford.edu/) and visualized by the “vioplot” package.

Establishment and evaluation of a pMRlncRNA risk signature

Using the R package “glmnet,” the least absolute shrinkage and selection operator (LASSO) regression analysis was implemented to further identify pMRlncRNAs that correlated with OS. Based on the variables obtained from LASSO and the regression coefficients, the risk score of each patient with melanoma was calculated using the following formula: risk score = sum of coefficients × pMRlncRNA expression. Then, we classified the patients into high- or low-risk groups based on the median risk score. Survival analysis was performed to determine whether a significant difference in OS existed between the 2 risk groups. The GSE65904 cohort was used for validation. Cox regression analyses were performed to evaluate whether our risk signature could serve as an independent prognostic factor for patients with melanoma. Using the R package “survivalROC,” the precision of the risk signature was validated by the area under the curve (AUC) of the receiver operating characteristic (ROC) curves. Gene set enrichment analysis (GSEA) was implemented to investigate the underlying biological pathways of the pMRlncRNAs with a simulation of 1,000 and false discovery rate (FDR) of 0.25.

Analyses of the risk score and clinicopathological characteristics

To further evaluate the stability of our pMRlncRNA risk signature, we performed survival analyses in groups with different clinicopathological characteristic (age, gender, TNM stage, and pathological stage). The relationship between risk score and clinicopathological factors was also investigated in TCGA cohort. We also explored the association between risk score and immune infiltration. Finally, we analyzed the expression level of PD-L1 in both groups.

Development of a nomogram for clinical use

Using the generalized linear signature regression algorithm, the risk score and other clinicopathological characteristics were used to develop a nomogram for prognosis. One-, 3-, and 5-year calibration curves were drawn to determine the predictive value of the nomogram by comparing the predicted survival events and the virtual observed outcomes.

Statistical analysis

R version 3.6.2 (R Core Team) was used to perform the statistical analysis. P<0.05 was statistically significant. Difference analysis between 2 clusters or groups was performed using Wilcoxon rank-sum test. Correlation analysis was performed by Pearson correlation analysis. Survival analysis was visualized by Kaplan-Meier curves and determined by log-rank test.

Results

Identification of prognostic pMRlncRNAs in melanoma

The workflow of this study is displayed in Figure 1. A total of 14,142 lncRNAs with expression data were extracted from TCGA database for subsequent analysis. The clinical characteristics of the included patients with melanoma are shown in Table 1. The expression matrices of 23 m6A-related genes were also generated. We performed Pearson correlation analysis to extract pMRlncRNAs with the criteria of P<0.001 and |R| >0.4. Finally, 234 pMRlncRNAs were obtained (Table S2). Fifteen prognostic pMRlncRNAs were identified by univariate Cox regression analysis (Figure 2A; P<0.05), all of which were protective.

Figure 1 Workflow of the study. TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic; PD-L1, programmed death-ligand 1.

Figure 2 Consensus clustering of pMRlncRNAs with prognostic value. (A) Forest plot shows the HR with a 95% confidence interval of 15 candidate prognosis-related pMRlncRNAs selected after univariate Cox regression analysis. (B) TCGA cohorts divided into two clusters based on the consensus clustering matrix (k=2). (C) Survival analysis demonstrates a better prognosis of patients from cluster 2. (D) Cluster 2 shows higher PD-L1 expression level than cluster 1. (E) Different distribution of clinicopathological characteristics in the two clusters. **, P<0.01, and ***, P<0.001. pMRlncRNAs, potential-N6-methyladenosine-related long noncoding RNAs; HR, hazard ratio; TCGA, The Cancer Genome Atlas; PD-L1, programmed cell death-Ligand 1.

Significant correlation of subtype clusters with clinical factors and survival

From k=2 to 9, k=2 was found to have ideal clustering stability based on the similarity displayed by pMRlncRNA expression levels (Figure 2B). Based on the expression levels, 447 patients with melanoma were divided into 2 clusters (cluster 1 =259; cluster 2 =188). Cluster 2 had a better OS and higher PD-L1 expression level than cluster 1 (Figure 2C,2D, both P<0.001). The expression levels of the pMRlncRNAs were higher in cluster 2 than in cluster 1 (Figure 2E). The clinical characteristics of the 2 clusters were also compared, and cluster 2 was significantly associated with a low T stage (Figure 2E, P<0.01).

The immune landscape in melanoma subtypes

Based on the results of CIBERSORT, the immune landscape of the 2 clusters was depicted, and the violin plot confirmed that there were significant differences in immune infiltration between the 2 clusters (Figure 3A). Cluster 1 had higher infiltration levels of M0 macrophages (P<0.05), M2 macrophages (P<0.05), and resting mast cells (P<0.01). However, cluster 2 had higher infiltration levels of memory B cells (P<0.001), plasma cells (P<0.01), resting memory CD4 T cells (P<0.05), activated memory CD4 T cells (P<0.05), T follicular helper cells (P<0.05), resting natural killer (NK) cells (P<0.05), and M1 macrophages (P<0.01). In addition, ESITIMATE results showed that cluster 2 had higher ESITIMATEScore (Figure 3B, P<0.001), ImmuneScore (Figure 3C, P<0.001), and StromalScore (Figure 3D, P<0.05) values than cluster 1.

Figure 3 Differences in immune infiltration in the two clusters. (A) Different infiltrating levels of 22 immune cell types in the 2 clusters. (B-D) Cluster 2 has higher ESTIMATEScore, ImmuneScore, and StromalScore than cluster 1.

Establishment and evaluation of the pMRlncRNAs risk signature

Using the R package “glmnet”, LASSO regression analysis was further implemented to identify pMRlncRNAs related to OS (Figure 4A,4B). Finally, 10 pMRlncRNAs with the highest prognostic value were identified (Table S3), and each patient’s risk score was obtained using the mentioned formula in TCGA cohort. Using the median as the cut-off value, all patients with melanoma were divided into high-risk and low-risk groups, based on their risk score. The patients are sorted by increasing risk score (Figure 4C). A scatter diagram showed the survival status of the patients (Figure 4D), while a heatmap demonstrated that the expression levels of the 10 pMRlncRNAs were relatively lower in the high-risk group (Figure 4E). Survival analysis confirmed that patients in the high-risk group had significantly worse OS than the low-risk group (Figure 4F, P<0.001). The above results were verified using the GSE65904 dataset (Figure 4G-4J). Finally, subgroup survival analyses were carried out based on the different clinical characteristics in TCGA cohort; we confirmed that our pMRlncRNA risk signature retained its disease prediction ability under different clinical conditions (Figure S1).

Figure 4 Development of a 10-pMRlncRNA risk signature for prognosis prediction. (A) Validation of the candidate pMRlncRNAs using LASSO regression analysis. (B) Explanation for LASSO coefficient profile plot of prognostic pMRlncRNAs. (C) Patients with melanoma patients sorted by growing risk score in TCGA cohort. (D) Survival status of patients with melanoma in TCGA cohort. (E) Heatmap showing the expression levels of 10 pMRlncRNA in the 2 risk groups from TCGA cohort. (F) Survival analysis performed for the 10-pMRlncRNA risk signature to assess overall survival in TCGA cohort. (G-J) Validation of the above results in the GSE65904 cohort. pMRlncRNAs, potential-N6-methyladenosine-related long noncoding RNAs; LASSO, least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas.

Independent prognostic analysis of the pMRlncRNAs risk signature

We further assessed the independence of the risk signature. Based on TCGA cohort, the results of univariate Cox analysis (HR: 2.933, 95% CI: 2.032–4.234, P<0.001; Figure 5A) and multivariate Cox analysis (HR: 2.342, 95% CI: 1.617–3.390, P<0.001; Figure 5B) showed that our risk signature was an independent prognostic factor for melanoma. The ROC curves also demonstrated that the risk score (AUC =0.720) had a preferable prognostic performance in predicting melanoma outcomes (Figure 5C). As shown in Figure 5D, 10 underlying pathways might be related to the development and outcome of melanoma.

Figure 5 Independent prognostic analysis and GSEA of the pMRlncRNAs risk signature. Univariate (A) and multivariate (B) Cox regression analyses to evaluate whether the risk score could be regarded as an independent prognostic factor of melanoma. (C) ROC curves reveal the AUC value of risk score and other clinicopathological characteristics. (D) GSEA of the 2 risk groups. GSEA, gene set enrichment analysis; ROC, receiver operating characteristic; AUC, area under the curve.

The risk score was associated with T stage, ImmunoScore, and cluster subtype in melanoma

We further explored the relationship between risk score and clinical characteristics. A heatmap demonstrated that the expression levels of the 10 pMRlncRNAs were upregulated in the low-risk group (Figure 6A). Significant differences in T stage (P<0.01), ImmuneScore (P<0.001), and cluster subtype (P<0.001) were observed between the 2 risk groups. Then, the relationship between ImmuneScore, T stage, and cluster subtype was analyzed in detail. The risk score increased with increase in T stage (P<0.001, Figure 6B). Cluster 1 had a higher risk score than cluster 2 (P<0.001, Figure 6C). The risk score of the high-ImmuneScore group was significantly lower than that of the low-ImmuneScore group (P<0.001, Figure 6D). The above results confirmed that the risk score was significantly correlated with T stage, Immunoscore, and cluster subtype in melanoma. Interestingly, the low-risk group also had significantly higher PD-L1 expression level than the high-risk group (P<0.001, Figure 6E).

Figure 6 The risk score correlates with clinicopathological factors in melanoma. (A) Heatmap showing significant differences in T stage, ImmuneScore, and cluster subtype between the 2 risk groups. The risk score in different (B) T stage, (C) cluster subtypes, and (D) ImmuneScore. (E) The low-risk group shows higher PD-L1 expression level than the high-risk group. **, P<0.01, and ***, P<0.001.

In addition, we also simultaneously explored the relationship between the risk score and immune infiltration, and found a negative correlation between the risk score and infiltration levels of M1 macrophages (Figure S2A, P<0.001), plasma cells (Figure S2B, P<0.001), activated memory CD4 T cells (Figure S2C, P<0.001), CD8 T cells (Figure S2D, P<0.001), T follicular helper cells (Figure S2E, P<0.001), and memory B cells (Figure S2F, P<0.01). Moreover, we found a positive correlation between the risk score and infiltration levels of M2 macrophages (Figure S2G, P<0.001), resting mast cells (Figure S2H, P<0.01), resting NK cells (Figure S2I, P<0.001), and M0 macrophages (Figure S2J, P<0.001).

Development of a nomogram for outcome prediction

A nomogram was built based on the risk score and other clinicopathological factors (Figure 7A). As shown by the 1-, 3-, and 5-year calibration curves (Figure 7B-7D), the nomogram was relatively accurate and stable in predicting melanoma prognosis.

Figure 7 Development of a nomogram combining risk score and clinicopathological features for prognosis prediction. (A) Nomogram combining risk score and clinicopathological features. (B-D) Calibration plots assessing the precision of the nomogram at predicting the 1-, 3-, and 5-year survival.

Discussion

Melanoma is the most dangerous form of skin cancer because of its high potential for distant metastasis and poor prognosis (24). Significant progress has been made in immunotherapy; however, its efficacy and safety vary for each patient. Thus, new biomarkers in melanoma are sorely needed for treatment decision-making and better prognosis. Previous studies have demonstrated that m6A methylation is the most common post-transcriptional RNA modification and affects the immune microenvironment (25). For instance, He et al. reported that the expression levels of m6A-related genes are significantly associated with prognosis and anti-tumor immune response in breast carcinoma (26). Zhou et al. developed an m6Ascore model to predict treatment response to ICIs and immune evasion in pancreatic carcinoma (27). Xu et al. found that the expression of the m6A eraser-related-genes, FTO and ALKBH5, is associated with gastric cancer prognosis (28). However, all the above studies were conducted at the genetic level. In addition, research focusing on the role of pMRlncRNAs in melanoma is still limited, and how pMRlncRNAs influence the immune microenvironment remains to be explored. In this study, we included TCGA and GSE65904 cohorts to determine the expression levels and prognostic significance of pMRlncRNAs in melanoma, and elucidate their effect on the tumor microenvironment.

A total of 15 protective pMRlncRNAs were identified through univariate Cox regression analysis. In addition, 2 clusters of patients with melanoma were identified through consensus clustering for pMRlncRNAs. Cluster 1 had a more advanced T stage and worse prognosis, whereas cluster 2 had higher PD-L1 expression level, and ESITIMATEScore, ImmuneScore, and StromalScore. The two clusters had different immune microenvironments. Cluster 2 had higher infiltration levels of memory B cells, plasma cells, CD4 T cells, NK cells, and M1 macrophages, which might explain the better OS in cluster 2. This step functionally elucidates the biological characteristics of the pMRlncRNA in melanoma. To reduce the overfitting and predict the clinical outcome of pMRlncRNA in melanoma patients precisely, LASSO regression analysis was performed. EBLN3P, SPRY4.AS1, WAC.AS1, AC093726.1, AC015911.3, AC083799.1, AC090559.1, C5orf56, AL133371.2, and AC092747.4 were finally included in the development of a pMRlncRNA risk signature, all of which were downregulated in the high-risk group.

The lncRNA EBLN3P is involved in the origin and progression of many diseases. Mathias et al. reported that EBLN3P is associated with immune response suppression and progression and related to good prognosis in the progression-free interval of breast cancer (29). Xu et al. reported that EBLN3P regulates UHMK1 expression by targeting miR-323a-3p/UHMK1 and accelerates colorectal cancer progression (30). Dai et al. demonstrated that EBLN3P accelerates the progression of osteosarcoma by regulating miR-224-5p/Rab10 and is a novel biomarker for the diagnosis and treatment of osteosarcoma (31). However, a literature search showed a limited number of studies reporting the biological mechanisms whereby the remaining pMRlncRNAs affect diseases. Therefore, our study may help identify vital pMRlncRNAs and provide new insights into their potential roles in melanoma tumorigenesis and progression.

Through survival analysis, we found that in TCGA cohort, high-risk patients had a significantly worse OS. Furthermore, this result was confirmed in the GSE65904 cohort.

Using Cox regression analyses, we confirmed that this pMRlncRNA risk signature is an independent prognostic factor for melanoma with high accuracy (AUC =0.720) that can be widely applied in subgroups with different clinical characteristics. To explore the underlying biological mechanism of the pMRlncRNAs, GSEA was implemented; 10 important Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to the tumorigenesis and progression of melanoma were listed. The JAK-STAT signaling pathway is an important therapeutic target for melanoma. Hu et al. reported that CXCL8 gene silencing may promote apoptosis of melanoma cells by inhibiting the JAK-STAT signaling pathway (32). Nguyen et al. also reported that mutations in the IFNγ-JAK-STAT pathway contributing to resistance to ICIs in melanoma enhance sensitivity to oncolytic viruses (33). The toll-like receptor signaling pathway also plays a vital role in melanoma. A previous study confirmed that toll-like receptor 4 signaling promotes the migration of melanoma, and compound A weakens toll-like receptor 4-mediated paclitaxel resistance in melanoma by suppressing interleukin-8 (34,35).

The relationship between the risk score and clinical characteristics was further investigated; we found that advanced T stages, cluster 1, and low ImmuneScore were associated with higher risk scores. In addition, patients with low-risk had higher PD-L1 expression level, which might explain the better OS in this group.

In addition, the risk score was combined with other clinicopathological characteristics, and a nomogram with calibration curves was developed, which demonstrated good predictive ability for predicting 1-, 3-, and 5-year survival, and the actual observed prognosis.

Our study has certain limitations. First, due to the lack of tumor samples, there are no biological experiments confirming our results. More functional studies on the 10 pMRlncRNAs should be conducted to further confirm the precision of the risk signature and elucidate the potential biological mechanisms. Second, multicenter and large-cohort studies are required to verify the risk signature before clinical use.

In conclusion, we identified 10 pMRlncRNAs related to melanoma OS, and developed a novel and independent pMRlncRNA risk signature for prognosis prediction. The advantage of this study is that this risk signature is based on high-throughput data from TCGA and GEO databases. These results may offer novel insights into the prognostic evaluation of melanoma and provide theoretical foundation for future studies on immune treatment in melanoma.

Acknowledgments

We thank Editage for the language editing and Dr. Nan Jiang and Dr. Zijian Yang for their helpful suggestions.

Funding: This study was supported by the National Natural Science Foundation of China (No. 81972559) and the Shanghai Shenkang Hospital Development Center Project (No. HDC2020CR2067B).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-21-4402/rc

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-21-4402/coif). JG reports that this study was supported by the National Natural Science Foundation of China (No. 81972559) and Shanghai Shenkang Hospital Development Center Project (No. HDC2020CR2067B). The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Xiao B, Liu L, Chen Z, et al. Identification of Epithelial-Mesenchymal Transition-Related Prognostic lncRNAs Biomarkers Associated With Melanoma Microenvironment. Front Cell Dev Biol 2021;9:679133. [Crossref] [PubMed]
2. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
3. Siegel RL, Miller KD, Fuchs HE, et al. Cancer Statistics, 2021. CA Cancer J Clin 2021;71:7-33. [Crossref] [PubMed]
4. Tsao H, Chin L, Garraway LA, et al. Melanoma: from mutations to medicine. Genes Dev 2012;26:1131-55. [Crossref] [PubMed]
5. Ho WJ, Erbe R, Danilova L, et al. Multi-omic profiling of lung and liver tumor microenvironments of metastatic pancreatic cancer reveals site-specific immune regulatory pathways. Genome Biol 2021;22:154. [Crossref] [PubMed]
6. Carter JM, Polley MC, Leon-Ferre RA, et al. Characteristics and Spatially Defined Immune (micro)landscapes of Early-stage PD-L1-positive Triple-negative Breast Cancer. Clin Cancer Res 2021;27:5628-37. [Crossref] [PubMed]
7. Yoneda S, Imagawa A, Hosokawa Y, et al. T-Lymphocyte Infiltration to Islets in the Pancreas of a Patient Who Developed Type 1 Diabetes After Administration of Immune Checkpoint Inhibitors. Diabetes Care 2019;42:e116-8. [Crossref] [PubMed]
8. Ribas A, Algazi A, Ascierto PA, et al. PD-L1 blockade in combination with inhibition of MAPK oncogenic signaling in patients with advanced melanoma. Nat Commun 2020;11:6262. [Crossref] [PubMed]
9. Wang Y, Li D, Lu J, et al. Long noncoding RNA TTN-AS1 facilitates tumorigenesis and metastasis by maintaining TTN expression in skin cutaneous melanoma. Cell Death Dis 2020;11:664. [Crossref] [PubMed]
10. Xie J, Zheng Y, Xu X, et al. Long Noncoding RNA CAR10 Contributes to Melanoma Progression By Suppressing miR-125b-5p to Induce RAB3D Expression. Onco Targets Ther 2020;13:6203-11. [Crossref] [PubMed]
11. Du K, Zhang L, Lee T, et al. m6A RNA Methylation Controls Neural Development and Is Involved in Human Diseases. Mol Neurobiol 2019;56:1596-606. [Crossref] [PubMed]
12. Dai D, Wang H, Zhu L, et al. N6-methyladenosine links RNA metabolism to cancer progression. Cell Death Dis 2018;9:124. [Crossref] [PubMed]
13. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol 2019;20:608-24. [Crossref] [PubMed]
14. Dahal U, Le K, Gupta M. RNA m6A methyltransferase METTL3 regulates invasiveness of melanoma cells by matrix metallopeptidase 2. Melanoma Res 2019;29:382-9. [Crossref] [PubMed]
15. Yang S, Wei J, Cui YH, et al. m6A mRNA demethylase FTO regulates melanoma tumorigenicity and response to anti-PD-1 blockade. Nat Commun 2019;10:2782. [Crossref] [PubMed]
16. Fu Y, Dominissini D, Rechavi G, et al. Gene expression regulation mediated through reversible m⁶A RNA methylation. Nat Rev Genet 2014;15:293-306. [Crossref] [PubMed]
17. Meyer KD, Jaffrey SR. Rethinking m6A Readers, Writers, and Erasers. Annu Rev Cell Dev Biol 2017;33:319-42. [Crossref] [PubMed]
18. Chelmicki T, Roger E, Teissandier A, et al. m6A RNA methylation regulates the fate of endogenous retroviruses. Nature 2021;591:312-6. [Crossref] [PubMed]
19. Edupuganti RR, Geiger S, Lindeboom RGH, et al. N6-methyladenosine (m6A) recruits and repels proteins to regulate mRNA homeostasis. Nat Struct Mol Biol 2017;24:870-8. [Crossref] [PubMed]
20. Arguello AE, DeLiberto AN, Kleiner RE. RNA Chemical Proteomics Reveals the N6-Methyladenosine (m6A)-Regulated Protein-RNA Interactome. J Am Chem Soc 2017;139:17249-52. [Crossref] [PubMed]
21. Shen S, Zhang R, Jiang Y, et al. Comprehensive analyses of m6A regulators and interactive coding and non-coding RNAs across 32 cancer types. Mol Cancer 2021;20:67. [Crossref] [PubMed]
22. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
23. Leek JT, Johnson WE, Parker HS, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012;28:882-3. [Crossref] [PubMed]
24. Yang Y, Long X, Li K, et al. Development and validation of an oxidative stress-associated prognostic risk model for melanoma. PeerJ 2021;9:e11258. [Crossref] [PubMed]
25. Ma S, Yan J, Barr T, et al. The RNA m6A reader YTHDF2 controls NK cell antitumor and antiviral immunity. J Exp Med 2021;218:e20210279. [Crossref] [PubMed]
26. He X, Tan L, Ni J, et al. Expression pattern of m6A regulators is significantly correlated with malignancy and antitumor immune response of breast cancer. Cancer Gene Ther 2021;28:188-96. [Crossref] [PubMed]
27. Zhou Z, Zhang J, Xu C, et al. An integrated model of N6-methyladenosine regulators to predict tumor aggressiveness and immune evasion in pancreatic cancer. EBioMedicine 2021;65:103271. [Crossref] [PubMed]
28. Xu X, Zhou E, Zheng J, et al. Prognostic and Predictive Value of m6A "Eraser" Related Gene Signature in Gastric Cancer. Front Oncol 2021;11:631803. [Crossref] [PubMed]
29. Mathias C, Muzzi JCD, Antunes BB, et al. Unraveling Immune-Related lncRNAs in Breast Cancer Molecular Subtypes. Front Oncol 2021;11:692170. [Crossref] [PubMed]
30. Xu XH, Song W, Li JH, et al. Long Non-coding RNA EBLN3P Regulates UHMK1 Expression by Sponging miR-323a-3p and Promotes Colorectal Cancer Progression. Front Med (Lausanne) 2021;8:651600. [Crossref] [PubMed]
31. Dai S, Li N, Zhou M, et al. LncRNA EBLN3P promotes the progression of osteosarcoma through modifying the miR-224-5p/Rab10 signaling axis. Sci Rep 2021;11:1992. [Crossref] [PubMed]
32. Hu X, Yuan L, Ma T. Mechanisms of JAK-STAT signaling pathway mediated by CXCL8 gene silencing on epithelial-mesenchymal transition of human cutaneous melanoma cells. Oncol Lett 2020;20:1973-81. [Crossref] [PubMed]
33. Nguyen TT, Ramsay L, Ahanfeshar-Adams M, et al. Mutations in the IFNγ-JAK-STAT Pathway Causing Resistance to Immune Checkpoint Inhibitors in Melanoma Increase Sensitivity to Oncolytic Virus Treatment. Clin Cancer Res 2021;27:3432-42. [Crossref] [PubMed]
34. Sootichote R, Thuwajit P, Singsuksawat E, et al. Compound A attenuates toll-like receptor 4-mediated paclitaxel resistance in breast cancer and melanoma through suppression of IL-8. BMC Cancer 2018;18:231. [Crossref] [PubMed]
35. Takazawa Y, Kiniwa Y, Ogawa E, et al. Toll-like receptor 4 signaling promotes the migration of human melanoma cells. Tohoku J Exp Med 2014;234:57-65. [Crossref] [PubMed]