Shichao Zhang1†Yu Xiong1†Chaochao Zheng2Jinhua Long2Houming Zhou1
Zhu Zeng1*Yan Ouyang2*
Fuzhou Tang2*- 1Key Laboratory of Infectious Immune and Antibody Engineering in Guizhou Province, School of Biology and Engineering, Guizhou Medical University, Guiyang, China
- 2Immune Cells and Antibody Engineering Research Center of Guizhou Province, School of Biology and Engineering, Guizhou Medical University, Guiyang, China
The “writers” of four types of adenosine (A)-related RNA modifications (N6-methyladenosine, N1-methyladenosine, alternative polyadenylation, as well as A-to-inosine RNA editing) are closely related to the tumorigenesis and progression of many cancer types, including skin cutaneous melanoma (SKCM). However, the potential roles of the crosstalk between these RNA modification “writers” in the tumor microenvironment (TME) remain unclear. The RNA modification patterns were identified using an unsupervised clustering method. Subsequently, based on differentially expressed genes responsible for the aforementioned RNA modification patterns, an RNA modification “writer” scoring model (W_Score) was constructed to quantify the RNA modification-associated subtypes in individual patients. Moreover, a correlation analysis for W_Score and the TME characteristics, clinical features, molecular subtypes, drug sensitivities, immune responses, and prognosis was performed. We identified three RNA modification patterns, corresponding to distinct tumor immune microenvironment characteristics and survival outcomes. Based on the W_Score score, which was extracted from the RNA modification-related signature genes, patients with SKCM were divided into high- and low-W_Score groups. The low-W_Score group was characterized by better survival outcomes and strengthened immunocyte infiltration. Further analysis showed that the low-W_Score group was positively associated with higher tumor mutation burden and PD-L1 expression. Of note, two immunotherapy cohorts demonstrated that patients with low W_Score exhibited long-term clinical benefits and an enhanced immune response. This study is the first to systematically analyze four types of A-related RNA modifications in SKCM, revealing that these “writers” essentially contribute to TME complexity and diversity. We quantitatively evaluated the RNA modification patterns in individual tumors, which could aid in developing personalized immunotherapy strategies for patients.
Introduction
Worldwide, there were 324,635 and 57,043 estimated new skin cutaneous melanoma (SKCM) cases and deaths, respectively, in 2020 (Sung et al., 2021). SKCM accounts for approximately 4% of all cases of skin cancer and is the most fatal subtype of skin cancer (Lin et al., 2021). Previous studies have demonstrated that the occurrence and development of SKCM are related to the accumulation of mutations in gene-modulating signaling pathways, including the Rb, p53, PI3K/AKT, and RAS/MAPK pathways (Dietrich et al., 2018; Chamcheu et al., 2019). Although the onset of SKCM may be partially attributed to somatic mutations, epigenetic changes in cancer-related genes are also associated with its etiology (de Unamuno Bustos et al., 2018). Epigenetic alterations mainly affect the functions and characteristics of genes by regulating gene transcription or translation (Bauer et al., 2016; Gu H. et al., 2020; Blanco et al., 2021). Furthermore, increasing evidence has revealed that RNA modification is mechanism that is indispensable for epigenetic regulation, which is involved in tumorigenesis and development of multiple cancers, including SKCM (Gu et al., 2015; Shen et al., 2019).
Over 170 types of RNA modifications have been identified (Wang et al., 2013; Barbieri and Kouzarides, 2020; Chen et al., 2021). Some of these modifications may interact to play key roles in many important biological processes, but synthetically analyzing all types of RNA modifications is difficult. Adenosine (A) is the most commonly modified nucleotide in RNA, with modifications that include N1-methyladenosine (m1A), N6-methyladenosine (m6A), and A-to-inosine (I) (Chen et al., 2021). Particularly, the m6A methylase negatively modulates the A-to-I RNA editing (Xiang et al., 2018). Thus, in this work, we focused on A-associated RNA modifications [m1A, m6A, A-to-I, and alternative polyadenylation (APA)] to explore the interaction by the “writers” that produced these modifications.
m6A is the most common RNA modification. Till date, this modification has been found in mRNAs, lncRNAs, miRNAs, circRNAs, rRNAs, tRNAs, and snRNAs (Zhao et al., 2017; Gu L. et al., 2020). Methyltransferases catalyze m6A methylation, and these “writers” include METTL14, METTL3, RBM15, RBM15B, VIRMA, WTAP, and ZC3H13 (Cao et al., 2016). In particular, m6A regulators with abnormal expression and genetic changes are involved in the pathogenesis and development of tumors, as well as in immune dysregulation (Li et al., 2017; Xiang et al., 2017; Zhang et al., 2017; Li et al., 2019; Gu C. et al., 2020). m1A is a methylation modification at the first nitrogen atom of the A base of the RNA molecule, which exists in tRNA, rRNA, and mRNA (Agris, 1996; Safra et al., 2017). The “writers” of the m1A modification include TRMT6, TRMT10C, TRMT61A, and TRMT61B (Dominissini et al., 2016; Wang Y. et al., 2019). Similar to the m6A modification, A-to-I RNA editing is also a prevailing RNA modification in the A base of mRNA (Di Giammartino et al., 2011; Abudayyeh et al., 2019). A-to-I RNA editing is mainly mediated by members of the A deaminase in the RNA (ADAR) family, including ADARB1, ADARB2, and ADAR, which bind to double-stranded RNA regions of protein-coding genes and non-coding sequences for the deamination of A to I (Huang et al., 2012; Han et al., 2014). APA is an important precursor-RNA processing mechanism widely present in all eukaryotes, which can regulate the length of the 3′-untranslated region (3′UTR) by cleaving mRNA at different sites, followed by the addition of poly(A) tails to the RNA 3′UTR (Elkon et al., 2013; Masamha et al., 2014). Several APA-related “writers” have been found to modulate the synthesis of poly (A) tails and the selection of variable poly (A) sites, including CFI, CLP1, CPSF1, CPSF2, CPSF3, CPSF4, CSTF1, CSTF2, CSTF3, NUDT21, PABPN1, and PCF11 (Batra et al., 2015).
The crosstalk among m1A, m6A, A-to-I, and APA can help reveal the mechanism underlying RNA modifications and the significance of post-transcriptional modifications. In SKCM, a complex regulatory network may be formed by the “writers” of the m1A, m6A, A-to-I, and APA modifications. Therefore, systematic analysis of the interaction between these “writers” could provide novel insights into the pathogenesis of SKCM and has potential clinical significance for tumor therapeutic target identification.
Current immunotherapies based on immune-checkpoint inhibitors (ICIs) have achieved astounding clinical efficacy in some patients. Unfortunately, immunotherapy is not effective in approximately 40–50% of patients with SKCM (Albittar et al., 2020). Tumor cells in the tumor microenvironment (TME) interact directly and indirectly with other components to induce hypoxia, chronic inflammation, and immunosuppression. Indeed, the abundance of tumor-infiltrating immune cells has been found to be associated with the prognosis of patients with SKCM (Ladányi et al., 2004). A high proportion of infiltrating CD8+ T cells appears to be a more effective immunotherapeutic response (Daud et al., 2016). Base on the characteristics of TME, tumors can be classified into immune-desert, immune-excluded, and immune-inflamed phenotypes, which present differences in the number of infiltrating immune cells and the immunotherapy response (Chen and Mellman, 2017; Zhang et al., 2020).
Recent studies have showed that the RNA modification “writers” are closely related to TME immune cell infiltration. Gao et al. revealed that m1A regulators, including “writers”, modulated the infiltration of immune cells in TME (Gao et al., 2021). Wang et al. found that m6A modification mediated by the RNA methyltransferase METTL3 promoted the functional activation of dendritic cells (DCs). Silencing METTL3 significantly downregulated the m6A modification levels, leading to a decreased antigen presentation ability of mature DCs (Wang H. et al., 2019). These studies were limited to one type of RNA modification, but multiple types of “writers” exert antitumor effects in a highly coordinated manner in cancer. Thus, systematic identification of the correlation between regulatory networks composed of multiple types of “writers” and the TME will be helpful in predicting immunotherapy responses and developing new immunotherapy strategies.
In this study, we comprehensively analyzed the relationship between four types of A-related RNA modification patterns and the characteristics of infiltrating immune cells by integrating transcriptomic data from 1526 SKCM samples from the Genotype-Tissue Expression (GTEx), Gene Expression Omnibus (GEO), and The Cancer Genome Atlas (TCGA) databases. Three distinct RNA modification patterns were identified, and the characteristics of the TME in these three patterns corresponded to immune-desert phenotype, immune-excluded phenotype, and immune-inflamed phenotype, respectively. We further constructed an RNA modification “writers” scoring model to quantify RNA modifications patterns in individual patients and to evaluate the patient response to targeted therapy and immunotherapy.
Methods
Data Acquisition of Skin Cutaneous Melanoma Samples
The workflow of this study is presented in Supplementary Figure S1A. Large-scale transcriptome data (including normal and tumor samples) was downloaded from the TCGA (https://tcga-data.nci.nih.gov/tcga/), GTEx (https://gtexportal.org/), and GEO (http://www.ncbi.nlm.nih.gov/geo) databases. Clinical information (including tumor stage, gender, age, and overall survival times) and somatic mutation and copy number variation (CNV) data from corresponding patients were retrieved from the University of California Santa Cruz genome browser (http://genome.ucsc.edu/). Patients with overall survival times <90 days and without survival information were excluded. Two GEO-SKCM cohorts (GSE78220 and GSE65904), one TCGA-SKCM cohort, and one GTEx-SK cohort were obtained for further analysis. The “limma with normalizeBetweenArrays” package was used to correct for batch effects of non-biotechnology deviation, and the basic information of these datasets is shown in Supplementary Table S1.
Unsupervised Clustering and Differential Expression Analysis of 26 “Writers”
Twenty-six RNA modification “writers” were obtained from the TCGA-SKCM cohort for identifying distinct RNA modification patterns. These 26 RNA modification “writers” included seven m6A regulators (METTL14, METTL3, RBM15, RBM15B, VIRMA, WTAP, and ZC3H13), four m1A regulators (TRMT6, TRMT10C, TRMT61A, and TRMT61B), three A-I regulators (ADARB1, ADARB2, and ADAR), and 12 APA regulators (CFI, CLP1, CPSF1, CPSF2, CPSF3, CPSF4, CSTF1, CSTF2, CSTF3, NUDT21, PABPN1, and PCF11). Based on the expression of 26 RNA modification “writers”, the tumor samples were classified (defined as “writer cluster”) via the unsupervised clustering analysis (Wilkerson and Hayes, 2010). For guaranteeing clustering stability, 1,000 iterations were performed with Spearman distance and pltem = 0.8 using the PAM algorithm (“ConsensuClusterPlus” package). The Wilcoxon signed-rank test was employed to analyze the differential expression of the 26 RNA modification “writers” between normal and tumor samples, and the expression correlation between these regulators was identified using Pearson’s correlation analysis.
Gene Set Variation Analysis and Functional Annotation
To examine the variations in RNA modification patterns in biological processes, gene set variation analysis (GSVA enrichment analysis; using the “GSVA” R package) was conducted (Hänzelmann et al., 2013). Gene sets with “c2.cp.kegg.v7.4.-symbols-gmt” and “c5.go.bp.v7.4.symbols.gmt” were retrieved from the Molecular Signatures Database followed by GSVA analysis with an adjusted p < 0.05. The functional annotation for 26 RNA modification “writers” was performed using the “clusterProfiler” R package (Yu et al., 2012).
Evaluation of Tumor-Infiltrating Immune Cells in the Tumor Microenvironment
Single-sample gene set enrichment analysis (ssGSEA) was employed to evaluate the relative proportion of each infiltrating immune cell type in the SKCM TME (Barbie et al., 2009; Charoentong et al., 2017), including adaptive immune cells (activated B cells, activated CD4+ T cells, activated CD8+ T cells, gamma delta T cells, immature B cells, regulatory T cells, T follicular helper cells, type 1 T helper cells, type 17 T helper cells, and type 2 T helper cells) and innate immune cells (activated DCs, CD56bright natural killer cells, CD56dim natural killer cells, eosinophils, immature DCs, macrophages, mast cells, MDSCs, monocytes, natural killer cells, natural killer T cells, neutrophils, and plasmacytoid DCs) (Supplementary Table S2). The enrichment scores were also calculated via ssGSEA.
Association Between Clusters and Other Biological Processes
Correlation analysis between categories and some associated biological processes was conducted. According to previous reports, gene sets including several gene-related biological processes were constructed, and these biological pathways included antigen processing and presentation, DNA damage repair, pan-fibroblast TGF-β response signature (pan-F-TBRS), angiogenesis, DNA replication, nucleotide excision repair, mismatch repair, WNT targets, cell cycle, CD8+ T effectors, antigen processing machinery, and immune-checkpoint genes (Şenbabaoğlu et al., 2016; Mariathasan et al., 2018).
Identification of Differentially Expressed Genes Between Different RNA Modification Phenotypes
To identify differentially expressed genes (DEGs) among the different RNA modification phenotypes (three distinct RNA modification patterns were identified based on the aforementioned unsupervised clustering analysis), an empirical Bayesian approach in the “limma” R package with an adjusted p < 0.05 was adopted (Ritchie et al., 2015). Furthermore, for functional annotation, we used the “clusterProfiler” R package to analyze DEGs with p < 0.05.
Constructing Gene Clusters
A univariate Cox regression analysis was performed to obtain prognosis-related DEGs, and unsupervised clustering analysis was then employed to classify these genes. In the PAM algorithm (using the “ConsensuClusterPlus” package) (Wilkerson and Hayes, 2010), 1,000 iterations were performed with Pearson distance and pltem = 0.8 to ensure the stability of the clusters.
Establishment of the W_Score Scoring System
Principal component analysis was used to construct the W_Score scoring model. First, we established a matrix consisting of expression levels of the prognosis-related DEGs of each patient. Then, principal components 1 and 2 were both selected act as signature scores. Finally, we calculated the W_Score for each cancer patient (similar to GGI) (Sotiriou et al., 2006):
where i denotes the expression of prognosis-related DEGs. Finally, according to the cutoff value (−0.646) determined using the “survminer” package, the patients with SKCM were divided into two categories (high W_Score and low W_Score).
Association Between W_Score and Immune-Related Signal Pathways
Correlation analysis between the W_Score and potential biological processes was conducted based on gene sets, including antigen processing and presentation, DNA damage repair, pan-F-TBRS, angiogenesis, DNA replication, nucleotide excision repair, mismatch repair, WNT targets, cell cycle, CD8+ T effectors, antigen processing machinery, immune-checkpoint, apoptosis regulation, ABL signaling, cytoskeleton, mitosis, ERK/MAPK signaling, RTK signaling, PI3K/MTOR_signaling, p53 pathway, protein stability and degradation, IGF1R signaling, genome integrity, JNK and p38 signaling, chromatin histone acetylation, EGFR signaling, metabolism, and WNT signaling pathways (Supplementary Table S6) (Iorio et al., 2016; Şenbabaoğlu et al., 2016; Mariathasan et al., 2018).
Evaluation of W_Score Correlations With Clinical Features and Drug Sensitivity
We evaluated the relationship between the clinicopathological characteristics, such as the Clark levels (I, II, III, IV, and V grades), transcriptomic classifications (immune, keratin, and MITF-low), and mutations subtypes (BRAF, RAS, and NF1) (Network, 2015), and the W_Score using the Wilcoxon signed-rank test. In addition, the data on drug response and targets/pathways were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) (http://www.cancerrxgene.org/downloads) (Supplementary Table S7) (Yang et al., 2013), and the Spearman’s correlation analysis was used to analyze the correlation between the W_Score and drug sensitivity. The “pRRophetic” package was used with an |Rs| > 0.2.
Analysis of W_Score and Immune-Checkpoint Genes
Gene expression profiles of the immune-checkpoint in patients with SKCM were retrieved from The Cancer Immunome Database (TCIA) (http://tcia.at/home) (Charoentong et al., 2017), including CTAL-4 positive and PD-1 negative (ips_ctla4_pos_pd1_neg), CTAL-4 positive and PD-1 positive (ips_ctla4_pos_pd1_pos), and CTAL-4 negative and PD-1 positive (ips_ctla4_neg_pd1_pos). Two tumor immunotherapeutic cohorts, including advanced urothelial cancer (treatment with atezolizumab, an anti-PD-L1 antibody) (Mariathasan et al., 2018), and metastatic melanoma (treatment with pembrolizumab, an anti-PD-1 antibody) (Hugo et al., 2016), were downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies and GEO (GSE78220), respectively.
Statistical Analysis
The differential expression analysis of 26 RNA modification “writers” between normal and tumor samples was conducted using the Wilcoxon signed-rank test. The correlation coefficients of the expression of 26 RNA modification “writers” were investigated using Pearson’s correlation analysis, while the correlation coefficients between tumor-infiltrating immune cells and 26 RNA modification “writers” were evaluated using the Wilcoxon signed-rank test. The hazard ratios (HRs) of RNA modification “writers” and DEGs were separately calculated via univariate Cox regression model. We used the Kaplan-Meier method with log-rank tests to obtain survival curves, and the reliability of the models was analyzed by the receiver operating characteristic (ROC) curve. A genome-wide CNV map with the 26 RNA modification “writers” was drawn using the R package “RCircos” (Mayakonda et al., 2018). All statistical analyses were two-sided, and p < 0.05 was considered statistically significant.
Results
Genetic Variation of Four Types of A-Associated RNA Modification “Writers” in SKCM
In this study, 26 A-associated RNA modification “writers”, including seven m6A modifications “writers”, four m1A modifications “writers”, three A-I modifications “writers”, and twelve APA modifications “writers”, were finally identified. We first evaluated the incidence of somatic mutations of the 26 RNA modification regulators in SKCM. Among the 467 samples, a mutation in at least one “writer” was observed in 130 patients (mutation frequency, 27.84%) (Figure 1A). Among the 26 “writers,” CFI exhibited the highest mutation frequency (6%), followed by VIRMA (5%) and ADARB2 (4%), whereas no mutations were found in TRMT61A and CPSF4 (Figure 1A). Notably, a significant co-occurrence between ZC3H13 and TRMT10C, ZC3H13 and CSTF1, ZC3H13 and ADAR, ADARB2 and RBM15, and VIRMA and PCF11 was observed (Supplementary Figure S2A). Furthermore, the patients were divided into two groups (the “writers” mutation group and the non-mutation group), and enrichment analysis via GSVA revealed that tumor hallmarks associated gene sets, including the KRAS signaling pathway, DNA repair pathway, and the PI3K/AKT/MTOR signaling pathway, were mainly enriched in “writers” mutation group (Supplementary Figure S2C). It has been demonstrated previously that the PI3K/AKT/MTOR signaling pathway is actively involved in the regulation of cancer cell proliferation, metastasis, survival, as well as the angiogenesis (Chamcheu et al., 2019). These findings suggested that mutations in A-associated RNA modification “writers” may cause functional changes thereby affecting SKCM progression.
FIGURE 1. The Genomic alterations and aberrant expression of RNA modification “writers” in SKCM. (A) The frequency of mutation of 26 RNA modification “writers” in 467 samples from the TCGA-SKCM cohort. (B) The location of copy number variation (CNV) of 26 RNA modification “writers” on 23 chromosomes. (C) The CNV frequency of 26 RNA modification “writers” in TCGA-SKCM cohort. (D) Principal component analysis based on the expression of 26 RNA modification “writers” to distinguish tumors (red dots) from normal tissues (blue dots) in the TCGA-SKCM and GTEx-SK cohorts. (E) The expression differences of 26 RNA modification “writers” between normal and tumor tissues. Tumor, red; Normal, blue. (*p < 0.05; **p < 0.01; ***p < 0.001).
We also investigated the CNVs of these “writers” and found that somatic copy number alterations were widespread in the 26 RNA modification regulator (Figure 1C). CPSF1 exhibited the highest frequency of CNV gains, followed by ADAR and CLP1, whereas CNV loss was observed in WTAP. Figure 1B showed the 26 RNA modification “writers” loci on a schematic of the genome. Interestingly, based on the expression levels of the 26 “writers”, the normal and tumor samples could be completely distinguished by principal component analysis (Figure 1D). We further discovered that the 26 RNA modification regulators have heterogeneous expression between normal and tumor tissues (Figure 1E), and the “writers” with amplified CNVs presented markedly higher expression in SKCM tissues, suggesting that CNV is one of the main factors that regulate “writers” expression (Supplementary Figure S3). The above analysis indicates that these RNA modification regulators with genomic alterations and expression imbalance play a potential role in the onset and development of SKCM.
Correlation Between RNA Modification “Writers” and TME Characteristics
The comprehensive landscape of the interactions of 26 RNA modification “writers” and their prognostic values in patients with SKCM was described in the regulator network (Figure 2A and Supplementary Figure S2B and Supplementary Table S3). We found that the expression of RNA regulators among the different types of RNA modification showed significant correlations. Therefore, crosstalk among these writers may play critical roles in the formation of distinct RNA modification patterns.
FIGURE 2. RNA modification patterns and corresponding biological characteristics of each pattern. (A) The correlations between 26 RNA modification “writers” expression in SKCM. The circle size represented the prognostic value of each regulator. (B) Survival analyses for the three RNA modification patterns based on 443 patients with SKCM from TCGA-SKCM cohort (Log-rank test; p = 0.007). (C–F) Comparison of GSVA and GO enrichment analysis (biological pathways) in distinct RNA modification patterns. (C,E), GSVA enrichment analysis; (D,F), GO enrichment analysis; (C,D), writer cluster A vs. writer cluster B; (E,F), writer cluster B vs. writer cluster C.
Based on the expression of the RNA modification “writers’’, unsupervised clustering analysis was then employed to stratify SKCM tumors with qualitatively different RNA modification patterns (Supplementary Figure S6A-D). A total of 3 clusters with distinct modification patterns were obtained, including 212 patients in cluster 1 (writer cluster A), 127 patients in cluster 2 (writer cluster B), and 103 patients in cluster 3 (writer cluster C). Next, based on the prognostic analysis for these clusters, patients in cluster_writer B were found to have a prominent survival advantage, whereas a poor prognosis was revealed for patients in writer cluster C (Figure 2B; log-rank test, p = 0.007). To further explore the biological significance underlying these distinct RNA modification patterns, GSVA enrichment analysis was conducted. The writer cluster A was significantly related to carcinogenic activation pathways (e.g., the TGF-β signaling pathway, the JAK-STAT signaling pathway, and the regulation of autophagy) (Figure 2C). The writer cluster B was mainly enriched in immune-associated signaling pathways (e.g., allograft rejection, chemokine signaling pathways, toll-like receptor signaling pathway, cytokine–cytokine receptor interactions, and T cell receptor signaling pathways) (Figure 2E), while writer cluster C was enriched in signaling pathways associated with metabolic reprogramming (e.g., tyrosine metabolism, glycine serine and threonine metabolism, and phenylalanine metabolism) (Figures 2C,E). These results were supported by gene ontology (GO; biological process) enrichment analysis (Figures 2D,F).
Because patients with writer cluster B presented a remarkable survival advantage and showed enrichment of immune-related pathways, we attempted to examine the functional roles of distinct RNA modification patterns in TME. The results showed that in the TCGA cohort (Supplementary Table S4), the writer cluster B exhibited increased fractions of tumor-infiltrating immune cells (e.g., mast cells, macrophages, MDSCs, plasmacytoid DCs, activated DCs, immature DCs, and natural killer cells) (Figure 3A), as well as a higher enrichment score for immune-related pathways (e.g., antigen processing and presentation) (Figure 3B), which was similar to the findings obtained for the GEO cohort (GSE65904; Supplementary Figure S4). Thus, writer cluster B could be characterized as an immune “hot” and activation phenotype, while writer cluster C was an immune “cold” and desert phenotype (Figures 3A,B). Previous studies have demonstrated that the prognosis for patients with SKCM is positively related to tumor-infiltrating immune cells, such as MDSCs, DCs, and natural killer cells, whereas it is negatively correlated with the infiltration of neutrophils and monocytes (Schmidt et al., 2007; Chevolet et al., 2015; Cursons et al., 2019), which is consistent with our results.
FIGURE 3. TME and transcriptome characteristics in distinct RNA modification patterns. (A) Evaluating the abundance of each tumor-infiltrating immune cell in three writer clusters using ssGSEA. (*p < 0.05; **p < 0.01; ***p < 0.001). (B) Differences in immune-associated pathways among three writer clusters (ANOVA analysis; *p < 0.05; **p < 0.01; ***p < 0.001). (C) Principal component analysis for the transcriptome profiles of patients in different writer clusters. (D) The association between RNA modification patterns (based on writer clusters) and clinical parameters, including age, gender, clinical stage, T-stage, N-stage, M-stage and survival status. (E) Functional annotation for differentially expressed genes (DEGs) (GO enrichment analysis, biological progress; BP).
Next, Spearman’s correlation analysis was conducted to investigate the specific associations between infiltrating immune cells and RNA modification “writers” (Supplementary Table S5A). We observed that WTAP (an m6A modification “writer”) demonstrated a positive correlation with a large number of TME infiltrating immune cells in SKCM. In particular, patients with elevated expression of WTAP presented remarkable more enrichment of TME DCs infiltration, including plasmacytoid DCs, activated DCs and immature DCs (Supplementary Figure S5B-C). Recently, Cao et al. proposed that RNA modulation “writers” can promote DC activation. DCs acting as antigen-presenting cells are a bridge between innate and adaptive immunities, and the activation of DCs relies on the high expression of adhesion factors, costimulatory factors, and MHC molecules (Qian and Cao, 2018). As expected, patients with a high expression of WTAP had an overall increase in the expression of adhesion factors, costimulatory factors, and MHC molecules (Supplementary Figure S5D). We also noted that the enhancement of immune-related pathways is accompanied by elevated expression of immune-checkpoint PD-1/L1 (CD279/CD274) in tumors with high expression of WTAP (Supplementary Figure S5D-E). Therefore, we evaluated the therapeutic effects of immune-checkpoint blocking antibodies between high and low WTAP expression patients. In anti-PD-1 immunotherapy cohort (GSE78220; Patients with metastatic melanoma treated with anti-PD-1 antibody immunotherapy), patients with a high expression of WTAP exhibited a survival benefit trend (Supplementary Figure S5F). Interestingly, WTAP was also associated with prolonged survival in GEO cohort (GSE65904; Supplementary Figure S5G). Taken together, WTAP may regulate the activation of DCs in the TME via RNA modification, participating in the antitumor immune response.
RNA Modification “Writer” Phenotype-Related DEGs
Despite patients with SKCM could be divided into three RNA modification phenotypes based on the expression of the RNA modification “writers” (Figures 3C,D and Supplementary Figure S6E), the genomic alterations and expression profile differences in these phenotypes are unclear. Therefore, we analyzed the changes in RNA modification-associated transcriptional expression across three RNA modification patterns. A total of 2281 RNA phenotype-associated DEGs were identified using an empirical Bayesian approach in the “limma” package (Supplementary Figure S6F). The “clusterProfiler” package was then used to perform GO analysis for these DEGs, and it showed that the biological processes were related with neutrophil-mediated immunity signaling pathway, neutrophil activation involved in immune response signaling pathway, regulation of the toll-like receptor signaling pathway, and the TRIF-dependent toll-like receptor signaling pathway (Figure 3E). 950 of these DEGs were significantly correlated with clinical outcomes, which were considered as RNA modification-associated signature genes. These results reconfirmed that RNA modification “writers” played an important role in the immune regulation in TME. Unsupervised clustering analysis was conducted to further verify this modulation mechanism, based on the gained 950 RNA phenotype-associated signature genes. This analysis classified patients into three genomic subtypes (Supplementary Figure S7A-D and Figure 4A), defined as gene cluster A, gene cluster B, and gene cluster C (183 cases, 83 cases, and 177 cases were clustered into gene cluster A, gene cluster B, and gene cluster C, respectively). We observed that the survival probability of patients was significantly different among these categories, and patients in gene cluster B presented the best prognosis compared with patients in gene cluster B and gene cluster C (with patients in gene cluster C showing the worst survival outcomes) (Figure 4C). Furthermore, a significant difference among the three categories was found for a variety of immune-related marker genes and enriched signaling pathways as well as the expression levels of 26 RNA modification “writers” (Figure 4B, and Supplementary Figure S7E-I). In addition, differences in prognosis and RNA “writer” expression between the three subgroups were further confirm by the analysis of patients with SKCM from GEO (GSE65904) (Supplementary Figure S8A-B). To our surprise, the characteristics of prognosis and the TME of gene cluster A, B, and C correspond to writer cluster A, B, and C, respectively.
FIGURE 4. Construction of the W_Score. (A) Unsupervised clustering of RNA modification phenotype-related differentially expressed genes (DEGs) in TCGA-SKCM cohort to classify patients into different genomic subtypes, defined as gene cluster A, B, and C. The writer clusters, gene clusters, age, gender, clinical stage, T-stage, N-stage, M-stage and survival status were used as patient annotations. (B) The differences in expression of 26 A-related RNA modification “writers” in three gene clusters (TCGA-SKCM cohort). (ANOVA analysis; *p < 0.05; **p < 0.01; ***p < 0.001). (C) The survival curves of three gene clusters based on 443 patients from TCGA-SKCM cohort (Log-rank test; p < 0.001). (D) The associations between W_Score and the known signaling pathways in TCGA-SKCM cohort. (E,F) Differences of the W_Score in writer clusters (E) and gene clusters (F), respectively, in TCGA-SKCM cohort (K-W test; p < 0.001). (G) Differences in immune-related pathways between high- and low- W_Score groups. (*p < 0.05; **p < 0.01; ***p < 0.001).
Generation of W_Score and Evaluation of Its Clinical Relevance
Considering the complexity and heterogeneity of RNA modification, we built a scoring model based on the RNA modification significant genes (950) to quantify the m6A modification pattern of individual patients with SKCM. We discovered that patients with writer cluster C exhibited a high W_Score, whereas those in writer cluster B di not. A similar association also observed between the W_Score and gene clusters (Figures 4E,F and Supplementary Table S5). Figure 5A showed the association among the writer cluster, Clark level, gene cluster, and W_Score. To evaluate the effect of the W_Score on TME, we compared the immune cell infiltration and scores of immune-related signaling pathway between the W_Score-low and -high groups. The results showed that patients with a low W_Score had a relatively high abundance of infiltrating immune cells in the TME (e.g., neutrophils, T follicular helper cells, MDSCs, activated DCs, mast cells, natural killer cells, macrophages, immature DCs, immature B cells, eosinophils, activated B cells, activated CD8+ T cells, activated CD4+ T cells, and regulatory T cells) (Supplementary Figure S9E) and a significant enhancement in immune activation pathways (e.g., CD8+ T effectors pathway, antigen processing machinery pathway, immune-checkpoints pathway, antigen processing and presentation pathway, and JNK and p38 signaling pathway) (Figures 4D,G). Therefore, tumors with low W_Score could be characterized by the immune “hot” phenotype, while high W_Score tumors could be closely linked to immune “cold” phenotype.
FIGURE 5. Clinical features associated with the W_Score. (A) Alluvial diagram exhibiting the relationship among the writer cluster, Clark level, gene cluster and W_Score. (B) Survival analyses for high- and low- W_Score groups in TCGA-SKCM cohort using Kaplan-Meier curves (Log-rank test; p < 0.001). (C–E) Distribution of W_Score in the different subtypes including gene mutations subtypes (C), immune subclasses (D), and the Clark levels (E). (F) Survival analyses for high- and low- TMB groups in TCGA-SKCM cohort (Log-rank test; p < 0.001). (G) Survival analyses for high TMB and high W_Score (H-TMB + H-W_Score), high TMB and low W_Score (H-TMB + L-W_Score), low TMB and high W_Score (L-TMB + H-W_Score), and low TMB and low W_Score (L-TMB + L-W_Score) patient groups in TCGA-SKCM cohort using Kaplan-Meier curves (Log-rank test; p < 0.001). (H–I) The waterfall plot of tumor somatic mutation constructed in low- (H) and high-W_Score (I) groups.
To further evaluate the clinical relevance of the W_Score, patients with SKCM were divided into high- and low-W_Score groups based on the cutoff value (0.0291) determined by the “survminer” package. Patients with a low W_Score had a survival advantage (Figure 5B; log-rank test, p < 0.001), and the 5-years survival rate was twice than patients in the high-W_Score group (46.1 vs 23.7%, respectively). Interestingly, multivariate Cox regression analysis (gender, age and T-stage, N-stage, M-stage, and clinical stage as covariates) confirmed that W_Score was an independent prognostic biomarker (HR = 1.022 (95% CI = 1.013–1.032), p < 0.001; Supplementary Figure S9A). The reliability of the W_Score was verified using 189 patients with SKCM from GEO (GSE65904). Consistent with the above findings, the W_Score-low group displayed a prominent survival benefit (Supplementary Figure S9F; log-rank test, p < 0.001), and W_Score could act as an independent prognostic indicator (HR = 1.082 (95% CI = 1.045–1.120), p < 0.001; Supplementary Figure S9B). These analyses indicate that W_Score can reflect the RNA modification patterns and be used to predict the clinical outcomes of patients with SKCM.
Subtypes distribution across different tumor grades and stages showed that patients diagnosed as advanced T-stage and higher Clark levels, as well as older ages had an elevated W_Score (Supplementary Figure S9C-D, and Figure 5E), implying the involvement of parameters comprising the W_Score in cancer progression. Based on the most markedly mutated genes, four subtypes were classified in the previous studies, including mutant BRAF, mutant RAS, mutant NF1, and mutant triple-wild-type (triple-WT, lacking hot-spot mutations in BRAF, RAS, or NF1). The BRAF subtype group presented the lowest W_Score (Figure 5C), suggestive of survival advantage, which was in accordance with the previous studies (Network, 2015). Moreover, SKCM patients in the keratin subclass were reported to exhibit worse prognosis than those in the immune subclass and the MITF-low subclass (Network, 2015). As expected, keratin subclass was significantly associated with a higher W_Score and worse prognosis (Figure 5D). These findings demonstrate the reliability of W_Score and augment the previous classification.
Increasing evidence has demonstrated that patients with high TMB have a favorable SKCM prognosis (Edwards et al., 2020; Yan et al., 2020). We also found that patients with a higher TMB displayed a survival advantage (Figure 5F). Therefore, we sought to identify the value of W_Score in evaluating the clinical outcome of patients with TMB. Patients with TMB were divided into four types, including high TMB and high W_Score (H-TMB + H-W_Score), high TMB and low W_Score (H-TMB + L-W_Score), low TMB and high W_Score (L-TMB + H-W_Score), and low TMB and low W_Score (L-TMB + L-W_Score). The patients in the H-TMB + L-W_Score group exhibited the best prognosis, whereas those in the L-TMB + H-W_Score group demonstrated the worst prognosis (Figure 5G). Furthermore, the TMB quantification analysis confirmed that a lower W_Score was remarkably related with a higher frequency of oncogene mutations, such as BRAF (56 vs. 45%), DNAH8 (32 vs. 24%), and DNAH10 (24 vs. 17%) (Figures 5H,I). These data enabled us more comprehensively to understand the effect of W_Score classification on genomic variation and prognosis of patients with SKCM, as well as to reveal the potential roles of RNA modification in the individual somatic mutations.
Potential Role of W_Score in Predicting Chemotherapeutic Drugs and Response to Immunotherapy With a PD-L1 Blocker
To further explore the effects of the W_Score on drug response, we evaluated the correlation between the W_Score and the response to drugs. A total of 48 drugs were identified from the GDSC database (Figure 6A). Among them, the sensitivity of 25 drugs was associated with the W_Score, including the farnesyltransferase inhibitor FTI-277 (Rs = −0.38, p = 6.72E-17), the WNT/β-catenin pathway inhibitor FH535 (Rs = −3.20, p = 9.07E-12), and the AKT and ERK inhibitors Sorafenib (Rs = −0.30, p = 2.93E-10), whereas 23 drugs showed the resistance related to the W_Score, including the DNA double-strand break and apoptosis activator Cisplatin (Rs = 0.38, p = 1.08E-16), the AKT inhibitor VIII (Rs = 0.28, p = 2.11E-09), and the cyclin-dependent kinase inhibitor CGP-60474 (Rs = 0.38, p = 1.01E-11). Then, the signaling pathway enrichment analysis based on drug-targeted genes revealed that the W_Score was linked to drug-associated signaling pathways, including the PI3K/MTOR signaling pathway, RTK signaling pathway, ERK/MAPK signaling pathway, p53 signaling pathway, and cell cycle signaling pathway (Figure 6B).
FIGURE 6. The W_Score predicting the drug sensitivity and immunotherapy efficacy. (A) The association between W_Score and drug sensitivity (Spearman analysis). (B) Signaling pathways targeted by drugs that are resistant (red) or sensitivity (blue) to the W_Score. (C,D) Differences in expressions of PD-L1 and CTAL-4 between high- and low- W_Score groups (Wilcoxon test; p < 0.0001). (E–G) The correlation between co-expression of PD-L1 and CTLA-4 and the W_Score (samples with SKCM from TCGA), including CTAL-4 positive and PD-1 negative (ips_ctla4_pos_pd1_neg) (E), CTAL-4 positive and PD-1 positive (ips_ctla4_pos_pd1_pos) (F), and CTAL-4 negative and PD-1 positive (ips_ctla4_neg_pd1_pos) (G). (H) Kaplan-Meier curve showing overall survival of SKCM patients between high- and low- W_Score groups in the anti-PD1 immunotherapy cohort (GSE78220 cohort; Log-rank test; p < 0.001). (I) The relationship of W_Score with clinical response to PD-1 blockade immunotherapy (GSE78220 cohort). (J) The proportion of patients with response to anti-PD-1 immunotherapy in high and low W_Score groups (GSE78220 cohort). PD, progressive disease; CR, complete response; PR, partial response. (K,L) The differences in the W_Score among distinct PD-1 blockade immunotherapy response groups (GSE78220 cohort). (M) The predictive value of W_Score in SKCM patients treated with anti-PD-1 immunotherapy (GSE78220; AUC, 0.863).
Immunotherapies based on CTLA-4/PD-1/PD-L1 inhibitors have achieved astounding clinical efficacy in malignant tumor therapy. Significant efforts have been made to identify biomarkers for predicting immunotherapy response, such as the expression of PD-L1 and CTLA-4 (Van Allen et al., 2015; Gibney et al., 2016; Miao et al., 2018). As shown in Figures 6C,D, the expression of PD-L1 and CTLA-4 obviously elevated in the low W_Score group. Subsequent analysis also revealed that a low W_Score was remarkably associated with high expression of PD-L1 and/or CTLA-4 (Figures 6E–G). Given that the W_Score was strongly related to the immune microenvironment, we further evaluated the power of the W_Score model in predicting patients’ responses to immunotherapy with PD-1/PD-L1 inhibitors, based on the anti-PD-1 cohort (GSE78220; metastatic melanoma) and anti-PD-L1 cohort (IMvigor210; advanced urothelial cancer). Patients with a lower W_Score displayed remarkably clinical benefits and a significantly prolonged survival in both cohorts (Figure 6H; p < 0.001, and Supplementary Figure S9G; p < 0.001). The better therapeutic outcomes and improved clinical response to anti-PD-1/L1 immunotherapy in W_Score-low patients were confirmed (Figures 6I–L, and Supplementary Figure S9H-I). Moreover, an AUC value with 0.870 indicated that the quantification of RNA modification patterns was a robust biomarker for evaluating patients’ prognosis and response to immunotherapy (Figure 6M). Taken together, these analyses imply that the W_score model could provide the selection of chemotherapeutic drugs and contribute to evaluating the response to anti-PD-1/L1 immunotherapy in SKCM.
Discussion
Mounting evidence has demonstrated that RNA modifications have an essential role in innate immunity, inflammation, and antitumor activity by interacting with multiple “writers” (Lin et al., 2020; Shulman and Stern-Ginossar, 2020; Chen et al., 2021; Chong et al., 2021). While current research has mainly focused on one type of RNA modification “writers”, such as m1A “writers” and m6A “writers” (Zhang et al., 2020; Gao et al., 2021), the roles of numberous types of RNA modification “writers” in tumor progression have not been fully recognized. Therefore, in this study, we explored global changes in four types of A-associated RNA modification “writers” (m6A, m1A, A-to-I, and APA) at the molecular level and the relationship between their interaction and the TME, so as to identify distinct RNA modification patterns in the tumor immune microenvironment. This would eventually reinforce our understanding of antitumor immune response and aid in the development of efficient immunotherapy strategies for patients with SKCM.
In this work, three distinct RNA modification patterns were identified based on 26 RNA modification “writers”. These three modification patterns showed remarkably different tumor immune microenvironment characteristics. The writer cluster B was characterized by the immune activation and abundant lymphocyte-infiltration, which correspond to the immune-inflamed phenotype. The writer cluster C was characterized by the immune-desert phenotype. The writer cluster A was characterized by the infiltration of immune cells together with the activation of TGF-β signaling pathway, JAK/STAT signaling pathway, and WNT signaling pathway, which correspond to the immune-excluded phenotype. The activated TGF-β pathways may suppress immune function by preventing lymphocytes from entering the tumor parenchyma, while specific molecular inhibitors that target TGF-β can restore the antitumor immunity by remodeling the tumor immune microenvironment (Fabregat et al., 2014; David et al., 2016; Mariathasan et al., 2018; Panagi et al., 2020). Studies have shown that the immune contexture of the TME influences the development and progression of SKCM, and can be used to predict immunotherapeutic response and prognosis (Galon and Bruni, 2019). The levels of tumor-infiltrating immune cells (e.g., natural killer cells, DCs, macrophages, and CD4 + /CD8 + T cells) are associated with the immune response that is generated (Topalian et al., 2016; Galon and Bruni, 2019; Zeng et al., 2020). Thus, these findings demonstrate the reliability of immunophenotypic classification based on distinct RNA modification patterns.
Furthermore, GO analysis and univariate Cox regression analysis suggested that prognosis-related DEGs could act as the signature genes of A-associated RNA modification. Of note, similar to the classification results of the RNA modification “writers”, the three gene categories based on A-associated RNA modification signature genes were constructed and were found to be remarkably related with different clinical outcomes and immunophenotypes. This confirmed that the significance of RNA modification “writers” in shaping the different landscape of TME. Because of the heterogeneity in RNA modification, we built a W_Score model to quantify the RNA modification patterns of individual patients. The RNA modification pattern characterized as the immune-inflamed phenotype displayed a lower W_Score, whereas those the modification patterns characterized as immune-desert and immune-excluded phenotypes presented a higher W_Score. We also found that the W_Score was significantly associated with clinical features of patients with SKCM and tumor subtypes. Noticeably, the W_Score was shown to be an independent prognostic biomarker for patients with SKCM, and this finding was verified in the GEO-SKCM cohort.
Because the benefits of survival and therapeutic responses of patients subjected to tumor immunotherapy remain limited to a small population, classifying patients to obtain better insight into the optimal use of tumor immunotherapy is an effective strategy (Harel et al., 2019; Albittar et al., 2020; Guo and Shen, 2020). Patients with SKCM were divided into high- and low-W_Score groups based on the cutoff value, and these two groups displayed extremely distinct tumor immune microenvironment characteristics. The low W_Score exhibited an immune “cold” phenotype characterized by a lack of immune-cell infiltration, indicative of a non-inflammatory TME (Kim and Chen, 2016). The high W_Score exhibited TME characteristics that were opposite to those of the immune “hot” phenotype, which was characterized by an enhanced immune-cell infiltration, indicative of an inflammatory tumor immune microenvironment (Turley et al., 2015; Chen and Mellman, 2017). These two groups could represent the distinct potential mechanisms modulating tumor immune escape, corresponding to different treatment strategies. The “hot” tumor (low W_Score group) exhibited a favorable immune-activated phenotype and may be responsive to ICIs (Wang et al., 2021). Notably, our data revealed that low W_Score showed a significantly high TMB. Analysis of genomic data indicated that a high TMB induced an increase in the number of neoantigens and boosted the immune response rate to ICI therapies (Liu et al., 2019; Hodi et al., 2021). In addition, patients with low W_Score had high enrichment scores in the immune-related signaling pathways. Roper et al. discovered that high expression of genes with antigen-processing machinery was beneficial for checkpoint immunotherapy (Roper et al., 2021). Ahn et al. reported that the PD-1 blockade promoted the activation of CD8+ T effectors and resulted in faster clearance of infection (Ahn et al., 2018). Furthermore, the clinical benefit of checkpoint immunotherapy is associated with high expression of immune-checkpoint proteins and suppression of angiogenesis (Castet et al., 2019; Zavareh et al., 2021). Of note, we confirmed that patients with low W_Score in two immunotherapy cohorts exhibited an enhanced immune response and long-term clinical benefits of immunotherapy. These analyses indicated that the W_Score could be regarded a powerful “predictor” for identifying patients sensitive to immunotherapy. In particular, the W_Score also has predictive values to evaluate the benefits of chemotherapy or targeted therapies. RNA modification patterns determined by the interaction of RNA modification “writers” were linked to immune phenotypes and affected the therapeutic effects of ICIs, which may provide clues on appropriate drug and immunotherapy strategies for SKCM.
In this study, we comprehensively analyzed RNA modification patterns based on four types of A-associated RNA modification “writers”, and revealed the association between these RNA modification patterns and tumor immune microenvironment characteristics. We also constructed a W_Score model and identified their clinical utility in targeted therapy and tumor immunotherapy. This paper has some limitations that further experimental verification are needed to validate these findings. This study thus highlighted the clinical significance of the crosstalk of multiple A-associated RNA modification. We believe that the findings of this study will contribute to developing personalized immunotherapy strategies for patients with SKCM.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.
Author Contributions
Conception and design: SZ and FT; Development of methodology: ZZ and YO. Acquisition of data: SZ, YX, ZZ, YO, and FT; Analysis and interpretation of data (e.g., statistical analysis, computational analysis): SZ, CZ, JL, HZ, and YX; Writing, review, and/or revision of the manuscript: SZ, ZZ, YO, and FT.
Funding
This study was supported by the National Natural Science Foundation of China (NSFC 31860244, 31760264, 32100442, and 31960139), the Science and Technology Foundation of Guizhou Province ((2021)172, (2020)1Z016, ZK (2021)025, (2020)1Y087, (2021)431, (2019)1275 and 19NSP002), Excellent Young Talents Plan of Guizhou Medical University (2020(105)), and the Science and Technology foundation of Guizhou Health Committee (gzwjkj 2019-1-037).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.821678/full#supplementary-material
Supplementary Figure S1 | Overview of study design.
Supplementary Figure S2 | Mutation characteristics of 26 A-related RNA modification “writers” in TCGA-SKCM cohort. (A) The mutation co-occurrence and exclusion analyses for 26 RNA regulators. Co-occurrence, blue; Exclusion, Khaki. (B) The prognostic analyses for 26 RNA regulators base on univariate Cox regression model. (C) GSVA enrichment analysis showed the differences between “writers” mutation and non-mutation tumors.
Supplementary Figure S3 | CNVs regulating the “writers” expression in TCGA-SKCM cohort. P value less than 0.05 indicated a statistical difference.
Supplementary Figure S4 | The characteristics of tumor microenvironment in three RNA modification patterns. Evaluating the abundance of each tumor-infiltrating immune cell in three writer clusters using ssGSEA (GSE65904). (*p < 0.05; **p < 0.01; ***p < 0.001).
Supplementary Figure S5 | Association between tumor immune microenvironment and 26 A-related RNA modification “writers”. (A) The correlation between the abundance of tumor-infiltrating immune cells and RNA regulators (TCGA-SKCM cohort; Spearman analysis; *p < 0.05; **p < 0.01; ***p < 0.001). (B) Differences of immuneScore between WTAP high expression and low expression groups (Wilcoxon test; p < 0.001). (C) Differences of tumor-infiltrating immune cells between WTAP high expression and low expression groups. (*p < 0.05; **p < 0.01; ***p < 0.001). (D) Differences in the expression of adhesion molecules, MHC molecules, and costimulatory molecules between WTAP high expression and low expression groups. (*p < 0.05; **p < 0.01; ***p < 0.001). (E) Differences in immune-activated pathways between WTAP high expression and low expression groups. (*p < 0.05; **p < 0.01; ***p < 0.001). (F,G) Kaplan-Meier curve showing overall survival of SKCM patients between high- and low- WTAP expression groups in the anti-PD-1 immunotherapy cohort (GSE78220 cohort) (F), and GEO-SKCM cohort (GSE65904) (G).
Supplementary Figure S6 | Unsupervised clustering of 26 A-related RNA modification “writers” in the TCGA-SKCM cohort. (A–D) Consensus matrices of the TCGA-SKCM cohort for k = 2–5. (E) The difference in expression of 26 RNA regulators among the writer cluster A, writer cluster B, and writer cluster C. (TCGA-SKCM cohort; *p < 0.05; **p < 0.01; ***p < 0.001). (F) 2281 differently expressed genes DEGs genes shown in venn diagram.
Supplementary Figure S7 | Immune-related molecular characteristics in distinct gene clusters. (A–D) Unsupervised clustering of prognosis-related differently expressed genes (DEGs) in the TCGA-SKCM cohort (k = 2–5). (E) Evaluating the abundance of each tumor-infiltrating immune cell in three gene clusters using ssGSEA. (*p < 0.05; **p < 0.01; ***p < 0.001). (F) Difference in immune-related pathways among three gene clusters. (*p < 0.05; **p < 0.01; ***p < 0.001). (G–I) Difference in the immune-activation (G), chemokines and cytokines (H), and immune-checkpoint (I) related gene expression among three gene clusters. (*p < 0.05; **p < 0.01; ***p < 0.001).
Supplementary Figure S8 | The expression of 26 “writers” in three gene clusters (GSE65904). (A) The survival curves of three gene clusters (GSE65904; Log-rank test; P < 0.05). (B) The differences in expression of 26 A-related RNA modification “writers” in three gene clusters. (GSE65904; ANOVA analysis; *p < 0.05; **p < 0.01; ***p < 0.001).
Supplementary Figure S9 | The correlation between W_Score and clinicopathological parameters. (A,B) Multivariate Cox regression analysis for W_Score in TCGA-SKCM cohort (A) and GEO-SKCM (GSE65904) cohort (B) shown by the forest plot. (C,D) Difference in W_Score among distinct clinical subgroups [age (C) and T-stage (D)] in TCGA-SKCM cohort. (*p < 0.05; **p <0.01; ***p < 0.001). (E) Evaluating the abundance of each tumor-infiltrating immune cell between high- and low- W_Score groups using ssGSEA. (F) Kaplan-Meier curve showing overall survival of SKCM patients between low- and high- W_Score groups (GSE65904) (Log-rank test; P < 0.001). (G) Kaplan-Meier curve showing overall survival of patients between low- and high- W_Score groups in the anti-PD-L1 immunotherapy cohort (IMvigor210 cohort; Log-rank test; P < 0.001). (H) The differences in the W_Score among distinct PD-L1 blockade immunotherapy response groups (IMvigor210 cohort). (I) The proportion of patients with response to PD-L1 blockade immunotherapy in high and low W_Score groups (IMvigor210 cohort). PD, progressive disease; CR, complete response; PR, partial response; SD, stable disease.
Abbreviations
3′UTR, 3′-untranslated region; APA, Alternative polyadenylation; A-related RNA modifications, adenosine-related RNA modifications; CNV, Copy number variations; DCs, Dendritic cells; DEGs, Differentially expressed genes; GDSC, Genomics of Drug Sensitivity in Cancer; GEO, Gene Expression Omnibus; GSVA, Gene Set Variation Analysis; GTEx, Genotype-Tissue Expression; HRs, hazard ratios; ICIs, immune-checkpoint inhibitors; ROC, receiver operating characteristic; SKCM, skin cutaneous melanoma; ssGSEA, Single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TCIA, The Cancer Immunome Database; TME, tumor microenvironment.
References
Abudayyeh, O. O., Gootenberg, J. S., Franklin, B., Koob, J., Kellner, M. J., Ladha, A., et al. (2019). A Cytosine Deaminase for Programmable Single-Base RNA Editing. Science 365 (6451), 382–386. doi:10.1126/science.aax7063
Agris, P. F. (1996). The Importance of Being Modified: Roles of Modified Nucleosides and Mg2+ in RNA Structure and Function. Prog. Nucleic Acid Res. Mol. Biol. 53, 79–129. doi:10.1016/s0079-6603(08)60143-9
Ahn, E., Araki, K., Hashimoto, M., Li, W., Riley, J. L., Cheung, J., et al. (2018). Role of PD-1 during Effector CD8 T Cell Differentiation. Proc. Natl. Acad. Sci. USA 115 (18), 4749–4754. doi:10.1073/pnas.1718217115
Albittar, A. A., Alhalabi, O., and Glitza Oliva, I. C. (2020). Immunotherapy for Melanoma. Adv. Exp. Med. Biol. 1244, 51–68. doi:10.1007/978-3-030-41008-7_3
Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA Interference Reveals that Oncogenic KRAS-Driven Cancers Require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460
Barbieri, I., and Kouzarides, T. (2020). Role of RNA Modifications in Cancer. Nat. Rev. Cancer 20 (6), 303–322. doi:10.1038/s41568-020-0253-2
Batra, R., Manchanda, M., and Swanson, M. S. (2015). Global Insights into Alternative Polyadenylation Regulation. RNA Biol. 12 (6), 597–602. doi:10.1080/15476286.2015.1040974
Bauer, T., Trump, S., Ishaque, N., Thürmann, L., Gu, L., Bauer, M., et al. (2016). Environment‐induced Epigenetic Reprogramming in Genomic Regulatory Elements in Smoking Mothers and Their Children. Mol. Syst. Biol. 12 (3), 861. doi:10.15252/msb.20156520
Blanco, M. A., Sykes, D. B., Gu, L., Wu, M., Petroni, R., Karnik, R., et al. (2021). Chromatin-state Barriers Enforce an Irreversible Mammalian Cell Fate Decision. Cel Rep. 37 (6), 109967. doi:10.1016/j.celrep.2021.109967
Cao, G., Li, H.-B., Yin, Z., and Flavell, R. A. (2016). Recent Advances in Dynamic M 6 A RNA Modification. Open Biol. 6 (4), 160003. doi:10.1098/rsob.160003
Castet, F., Garcia-Mulero, S., Sanz-Pamplona, R., Cuellar, A., Casanovas, O., Caminal, J., et al. (2019). Uveal Melanoma, Angiogenesis and Immunotherapy, Is There Any hope? Cancers 11 (6), 834. doi:10.3390/cancers11060834
Chamcheu, J., Roy, T., Uddin, M., Banang-Mbeumi, S., Chamcheu, R.-C., Walker, A., et al. (2019). Role and Therapeutic Targeting of the PI3K/Akt/mTOR Signaling Pathway in Skin Cancer: A Review of Current Status and Future Trends on Natural and Synthetic Agents Therapy. Cells 8 (8), 803. doi:10.3390/cells8080803
Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019
Chen, D. S., and Mellman, I. (2017). Elements of Cancer Immunity and the Cancer-Immune Set point. Nature 541 (7637), 321–330. doi:10.1038/nature21349
Chen, H., Yao, J., Bao, R., Dong, Y., Zhang, T., Du, Y., et al. (2021). Cross-talk of Four Types of RNA Modification Writers Defines Tumor Microenvironment and Pharmacogenomic Landscape in Colorectal Cancer. Mol. Cancer 20 (1), 29. doi:10.1186/s12943-021-01322-w
Chevolet, I., Speeckaert, R., Schreuer, M., Neyns, B., Krysko, O., Bachert, C., et al. (2015). Clinical Significance of Plasmacytoid Dendritic Cells and Myeloid-Derived Suppressor Cells in Melanoma. J. Transl Med. 13, 9. doi:10.1186/s12967-014-0376-x
Chong, W., Shang, L., Liu, J., Fang, Z., Du, F., Wu, H., et al. (2021). m6A Regulator-Based Methylation Modification Patterns Characterized by Distinct Tumor Microenvironment Immune Profiles in colon Cancer. Theranostics 11 (5), 2201–2217. doi:10.7150/thno.52717
Cursons, J., Souza-Fonseca-Guimaraes, F., Foroutan, M., Anderson, A., Hollande, F., Hediyeh-Zadeh, S., et al. (2019). A Gene Signature Predicting Natural Killer Cell Infiltration and Improved Survival in Melanoma Patients. Cancer Immunol. Res. 7 (7), 1162–1174. doi:10.1158/2326-6066.cir-18-0500
Daud, A. I., Loo, K., Pauli, M. L., Sanchez-Rodriguez, R., Sandoval, P. M., Taravati, K., et al. (2016). Tumor Immune Profiling Predicts Response to Anti-PD-1 Therapy in Human Melanoma. J. Clin. Invest. 126 (9), 3447–3452. doi:10.1172/jci87324
David, C. J., Huang, Y.-H., Chen, M., Su, J., Zou, Y., Bardeesy, N., et al. (2016). TGF-β Tumor Suppression through a Lethal EMT. Cell 164 (5), 1015–1030. doi:10.1016/j.cell.2016.01.009
de Unamuno Bustos, B., Murria Estal, R., Pérez Simó, G., Simarro Farinos, J., Pujol Marco, C., Navarro Mira, M., et al. (2018). Aberrant DNA Methylation Is Associated with Aggressive Clinicopathological Features and Poor Survival in Cutaneous Melanoma. Br. J. Dermatol. 179 (2), 394–404. doi:10.1111/bjd.16254
Dietrich, P., Kuphal, S., Spruss, T., Hellerbrand, C., and Bosserhoff, A. K. (2018). Wild-type KRAS Is a Novel Therapeutic Target for Melanoma Contributing to Primary and Acquired Resistance to BRAF Inhibition. Oncogene 37 (7), 897–911. doi:10.1038/onc.2017.391
Di Giammartino, D. C., Nishida, K., and Manley, J. L. (2011). Mechanisms and Consequences of Alternative Polyadenylation. Mol. Cel 43 (6), 853–866. doi:10.1016/j.molcel.2011.08.017
Dominissini, D., Nachtergaele, S., Moshitch-Moshkovitz, S., Peer, E., Kol, N., Ben-Haim, M. S., et al. (2016). The Dynamic N1-Methyladenosine Methylome in Eukaryotic Messenger RNA. Nature 530 (7591), 441–446. doi:10.1038/nature16998
Edwards, J., Ferguson, P. M., Lo, S. N., Pires da Silva, I., Colebatch, A. J., Lee, H., et al. (2020). Tumor Mutation burden and Structural Chromosomal Aberrations Are Not Associated with T-Cell Density or Patient Survival in Acral, Mucosal, and Cutaneous Melanomas. Cancer Immunol. Res. 8 (11), 1346–1353. doi:10.1158/2326-6066.cir-19-0835
Elkon, R., Ugalde, A. P., and Agami, R. (2013). Alternative Cleavage and Polyadenylation: Extent, Regulation and Function. Nat. Rev. Genet. 14 (7), 496–506. doi:10.1038/nrg3482
Fabregat, I., Fernando, J., Mainez, J., and Sancho, P. (2014). TGF-beta Signaling in Cancer Treatment. Curr. Pharm. Des. 20 (17), 2934–2947. doi:10.2174/13816128113199990591
Galon, J., and Bruni, D. (2019). Approaches to Treat Immune Hot, Altered and Cold Tumours with Combination Immunotherapies. Nat. Rev. Drug Discov. 18 (3), 197–218. doi:10.1038/s41573-018-0007-y
Gao, Y., Wang, H., Li, H., Ye, X., Xia, Y., Yuan, S., et al. (2021). Integrated Analyses of m1A Regulator-Mediated Modification Patterns in Tumor Microenvironment-Infiltrating Immune Cells in colon Cancer. Oncoimmunology 10 (1), 1936758. doi:10.1080/2162402x.2021.1936758
Gibney, G. T., Weiner, L. M., and Atkins, M. B. (2016). Predictive Biomarkers for Checkpoint Inhibitor-Based Immunotherapy. Lancet Oncol. 17 (12), e542–e551. doi:10.1016/S1470-2045(16)30406-5
Gu, C., Shi, X., Dai, C., Shen, F., Rocco, G., Chen, J., et al. (2020). RNA m6A Modification in Cancers: Molecular Mechanisms and Potential Clinical Applications. The Innovation 1 (3), 100066. doi:10.1016/j.xinn.2020.100066
Gu, H., Guo, Y., Gu, L., Wei, A., Xie, S., Ye, Z., et al. (2020). Deep Learning for Identifying Corneal Diseases from Ocular Surface Slit-Lamp Photographs. Sci. Rep. 10 (1), 17851. doi:10.1038/s41598-020-75027-3
Gu, L., Frommel, S. C., Frommel, S. C., Oakes, C. C., Simon, R., Grupp, K., et al. (2015). BAZ2A (TIP5) Is Involved in Epigenetic Alterations in Prostate Cancer and its Overexpression Predicts Disease Recurrence. Nat. Genet. 47 (1), 22–30. doi:10.1038/ng.3165
Gu, L., Wang, L., Chen, H., Hong, J., Shen, Z., Dhall, A., et al. (2020). CG14906 (Mettl4) Mediates m6A Methylation of U2 snRNA in Drosophila. Cell Discov 6, 44. doi:10.1038/s41421-020-0178-7
Guo, X., and Shen, W. (2020). Latest Evidence on Immunotherapy for Cholangiocarcinoma. Oncol. Lett. 20 (6), 381. doi:10.3892/ol.2020.12244
Han, S.-W., Kim, H.-P., Shin, J.-Y., Jeong, E.-G., Lee, W.-C., Kim, K. Y., et al. (2014). RNA Editing in RHOQ Promotes Invasion Potential in Colorectal Cancer. J. Exp. Med. 211 (4), 613–621. doi:10.1084/jem.20132209
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Harel, M., Ortenberg, R., Varanasi, S. K., Mangalhara, K. C., Mardamshina, M., Markovits, E., et al. (2019). Proteomics of Melanoma Response to Immunotherapy Reveals Mitochondrial Dependence. Cell 179 (1), 236–250. doi:10.1016/j.cell.2019.08.012
Hodi, F. S., Wolchok, J. D., Schadendorf, D., Larkin, J., Long, G. V., Qian, X., et al. (2021). TMB and Inflammatory Gene Expression Associated with Clinical Outcomes Following Immunotherapy in Advanced Melanoma. Cancer Immunol. Res. 9 (10), 1202–1213. doi:10.1158/2326-6066.cir-20-0983
Huang, H., Tan, B. Z., Shen, Y., Tao, J., Jiang, F., Sung, Y. Y., et al. (2012). RNA Editing of the IQ Domain in Cav1.3 Channels Modulates Their Ca2+-dependent Inactivation. Neuron 73 (2), 304–316. doi:10.1016/j.neuron.2011.11.022
Hugo, W., Zaretsky, J. M., Sun, L., Song, C., Moreno, B. H., Hu-Lieskovan, S., et al. (2016). Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 165 (1), 35–44. doi:10.1016/j.cell.2016.02.065
Iorio, F., Knijnenburg, T. A., Vis, D. J., Bignell, G. R., Menden, M. P., Schubert, M., et al. (2016). A Landscape of Pharmacogenomic Interactions in Cancer. Cell 166 (3), 740–754. doi:10.1016/j.cell.2016.06.017
Kim, J. M., and Chen, D. S. (2016). Immune Escape to PD-L1/pd-1 Blockade: Seven Steps to success (Or Failure). Ann. Oncol. 27 (8), 1492–1504. doi:10.1093/annonc/mdw217
Ladányi, A., Somlai, B., Gilde, K., Fejös, Z., Gaudi, I., and Tímár, J. (2004). T-cell Activation Marker Expression on Tumor-Infiltrating Lymphocytes as Prognostic Factor in Cutaneous Malignant Melanoma. Clin. Cancer Res. 10 (2), 521–530. doi:10.1158/1078-0432.ccr-1161-03
Li, H.-B., Tong, J., Zhu, S., Batista, P. J., Duffy, E. E., Zhao, J., et al. (2017). m6A mRNA Methylation Controls T Cell Homeostasis by Targeting the IL-7/STAT5/SOCS Pathways. Nature 548 (7667), 338–342. doi:10.1038/nature23450
Li, T., Hu, P.-S., Zuo, Z., Lin, J.-F., Li, X., Wu, Q.-N., et al. (2019). METTL3 Facilitates Tumor Progression via an m6A-igf2bp2-dependent Mechanism in Colorectal Carcinoma. Mol. Cancer 18 (1), 112. doi:10.1186/s12943-019-1038-7
Lin, Y., Wang, S., Liu, S., Lv, S., Wang, H., and Li, F. (2021). Identification and Verification of Molecular Subtypes with Enhanced Immune Infiltration Based on m6A Regulators in Cutaneous Melanoma. Biomed. Res. Int. 2021, 2769689. doi:10.1155/2021/2769689
Lin, Z., Niu, Y., Wan, A., Chen, D., Liang, H., Chen, X., et al. (2020). RNA M6 A Methylation Regulates Sorafenib Resistance in Liver Cancer through FOXO3-Mediated Autophagy. EMBO J. 39 (12), e103181. doi:10.15252/embj.2019103181
Liu, L., Bai, X., Wang, J., Tang, X.-R., Wu, D.-H., Du, S.-S., et al. (2019). Combination of TMB and CNA Stratifies Prognostic and Predictive Responses to Immunotherapy across Metastatic Cancer. Clin. Cancer Res. 25 (24), 7413–7423. doi:10.1158/1078-0432.ccr-19-0558
Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501
Masamha, C. P., Xia, Z., Yang, J., Albrecht, T. R., Li, M., Shyu, A.-B., et al. (2014). CFIm25 Links Alternative Polyadenylation to Glioblastoma Tumour Suppression. Nature 510 (7505), 412–416. doi:10.1038/nature13261
Mayakonda, A., Lin, D.-C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118
Miao, D., Margolis, C. A., Gao, W., Voss, M. H., Li, W., Martini, D. J., et al. (2018). Genomic Correlates of Response to Immune Checkpoint Therapies in clear Cell Renal Cell Carcinoma. Science 359 (6377), 801–806. doi:10.1126/science.aan5951
Network, C. G. A. (2015). Genomic Classification of Cutaneous Melanoma. Cell 161 (7), 1681–1696. doi:10.1016/j.cell.2015.05.044
Panagi, M., Voutouri, C., Mpekris, F., Papageorgis, P., Martin, M. R., Martin, J. D., et al. (2020). TGF-β Inhibition Combined with Cytotoxic Nanomedicine Normalizes Triple Negative Breast Cancer Microenvironment towards Anti-tumor Immunity. Theranostics 10 (4), 1910–1922. doi:10.7150/thno.36936
Qian, C., and Cao, X. (2018). Dendritic Cells in the Regulation of Immunity and Inflammation. Semin. Immunol. 35, 3–11. doi:10.1016/j.smim.2017.12.002
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Roper, N., Velez, M. J., Chiappori, A., Kim, Y. S., Wei, J. S., Sindiri, S., et al. (2021). Notch Signaling and Efficacy of PD-1/pd-L1 Blockade in Relapsed Small Cell Lung Cancer. Nat. Commun. 12 (1), 3880. doi:10.1038/s41467-021-24164-y
Safra, M., Sas-Chen, A., Nir, R., Winkler, R., Nachshon, A., Bar-Yaacov, D., et al. (2017). The m1A Landscape on Cytosolic and Mitochondrial mRNA at Single-Base Resolution. Nature 551 (7679), 251–255. doi:10.1038/nature24456
Schmidt, H., Suciu, S., Punt, C. J. A., Gore, M., Kruit, W., Patel, P., et al. (2007). Pretreatment Levels of Peripheral Neutrophils and Leukocytes as Independent Predictors of Overall Survival in Patients with American Joint Committee on Cancer Stage IV Melanoma: Results of the EORTC 18951 Biochemotherapy Trial. J. Clin. Oncol. 25 (12), 1562–1569. doi:10.1200/jco.2006.09.0274
Şenbabaoğlu, Y., Gejman, R. S., Winer, A. G., Liu, M., Van Allen, E. M., de Velasco, G., et al. (2016). Tumor Immune Microenvironment Characterization in clear Cell Renal Cell Carcinoma Identifies Prognostic and Immunotherapeutically Relevant Messenger RNA Signatures. Genome Biol. 17 (1), 231. doi:10.1186/s13059-016-1092-z
Shen, S., Faouzi, S., Bastide, A., Martineau, S., Malka-Mahieu, H., Fu, Y., et al. (2019). An Epitranscriptomic Mechanism Underlies Selective mRNA Translation Remodelling in Melanoma Persister Cells. Nat. Commun. 10 (1), 5713. doi:10.1038/s41467-019-13360-6
Shulman, Z., and Stern-Ginossar, N. (2020). The RNA Modification N6-Methyladenosine as a Novel Regulator of the Immune System. Nat. Immunol. 21 (5), 501–512. doi:10.1038/s41590-020-0650-4
Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade to Improve Prognosis. J. Natl. Cancer Inst. 98 (4), 262–272. doi:10.1093/jnci/djj052
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A. Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Topalian, S. L., Taube, J. M., Anders, R. A., and Pardoll, D. M. (2016). Mechanism-driven Biomarkers to Guide Immune Checkpoint Blockade in Cancer Therapy. Nat. Rev. Cancer 16 (5), 275–287. doi:10.1038/nrc.2016.36
Turley, S. J., Cremasco, V., and Astarita, J. L. (2015). Immunological Hallmarks of Stromal Cells in the Tumour Microenvironment. Nat. Rev. Immunol. 15 (11), 669–682. doi:10.1038/nri3902
Van Allen, E. M., Miao, D., Schilling, B., Shukla, S. A., Blank, C., Zimmer, L., et al. (2015). Genomic Correlates of Response to CTLA-4 Blockade in Metastatic Melanoma. Science 350 (6257), 207–211. doi:10.1126/science.aad0095
Wang, H., Hu, X., Huang, M., Liu, J., Gu, Y., Ma, L., et al. (2019). Mettl3-mediated mRNA m6A Methylation Promotes Dendritic Cell Activation. Nat. Commun. 10 (1), 1898. doi:10.1038/s41467-019-09903-6
Wang, Q., Gu, L., Adey, A., Radlwimmer, B., Wang, W., Hovestadt, V., et al. (2013). Tagmentation-based Whole-Genome Bisulfite Sequencing. Nat. Protoc. 8 (10), 2022–2032. doi:10.1038/nprot.2013.118
Wang, Y., Huang, Q., Deng, T., Li, B.-H., and Ren, X.-Q. (2019). Clinical Significance of TRMT6 in Hepatocellular Carcinoma: A Bioinformatics-Based Study. Med. Sci. Monit. 25, 3894–3901. doi:10.12659/msm.913556
Wang, Z., Little, N., Chen, J., Lambesis, K. T., Le, K. T., Han, W., et al. (2021). Immunogenic Camptothesome Nanovesicles Comprising Sphingomyelin-Derived Camptothecin Bilayers for Safe and Synergistic Cancer Immunochemotherapy. Nat. Nanotechnol 16 (10), 1130–1140. doi:10.1038/s41565-021-00950-z
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170
Xiang, J.-F., Yang, Q., Liu, C.-X., Wu, M., Chen, L.-L., and Yang, L. (2018). N6-methyladenosines Modulate A-To-I RNA Editing. Mol. Cel 69 (1), 126–135. doi:10.1016/j.molcel.2017.12.006
Xiang, Y., Laurent, B., Hsu, C.-H., Nachtergaele, S., Lu, Z., Sheng, W., et al. (2017). RNA m6A Methylation Regulates the Ultraviolet-Induced DNA Damage Response. Nature 543 (7646), 573–576. doi:10.1038/nature21671
Yan, J., Wu, X., Yu, J., Zhu, Y., and Cang, S. (2020). Prognostic Role of Tumor Mutation burden Combined with Immune Infiltrates in Skin Cutaneous Melanoma Based on Multi-Omics Analysis. Front. Oncol. 10, 570654–654. doi:10.3389/fonc.2020.570654
Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of Drug Sensitivity in Cancer (GDSC): a Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res. 41 (Database issue), D955–D961. doi:10.1093/nar/gks1111
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). ClusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS. 16 (5), 284–287. doi:10.1089/omi.2011.0118
Zavareh, R. B., Spangenberg, S. H., Woods, A., Martínez-Peña, F., and Lairson, L. L. (2021). HSP90 Inhibition Enhances Cancer Immunotherapy by Modulating the Surface Expression of Multiple Immune Checkpoint Proteins. Cel Chem. Biol. 28 (2), 158–168. doi:10.1016/j.chembiol.2020.10.005
Zeng, D., Ye, Z., Wu, J., Zhou, R., Fan, X., Wang, G., et al. (2020). Macrophage Correlates with Immunophenotype and Predicts Anti-PD-L1 Response of Urothelial Cancer. Theranostics 10 (15), 7002–7014. doi:10.7150/thno.46176
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0
Zhang, C., Chen, Y., Sun, B., Wang, L., Yang, Y., Ma, D., et al. (2017). m6A Modulates Haematopoietic Stem and Progenitor Cell Specification. Nature 549 (7671), 273–276. doi:10.1038/nature23883
Keywords: RNA modification “writer”, skin cutaneous melanoma, tumor microenvironment, W_Score, immunotherapy
Citation: Zhang S, Xiong Y, Zheng C, Long J, Zhou H, Zeng Z, Ouyang Y and Tang F (2022) Crosstalk Between Four Types of RNA Modification Writers Characterizes the Tumor Immune Microenvironment Infiltration Patterns in Skin Cutaneous Melanoma. Front. Cell Dev. Biol. 10:821678. doi: 10.3389/fcell.2022.821678
Received: 24 November 2021; Accepted: 10 January 2022;
Published: 26 January 2022.
Edited by:
Tao Huang, Shanghai Institute of Nutrition and Health (CAS), ChinaReviewed by:
Lei Gu, Max Planck Institute for Heart and Lung Research, GermanyChang Gu, Tongji University, China
Copyright © 2022 Zhang, Xiong, Zheng, Long, Zhou, Zeng, Ouyang and Tang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fuzhou Tang, tangfuzhou@163.com; Yan Ouyang, ouyangyan@gmc.edu.cn; Zhu Zeng, 724914301@qq.com
†These authors have contributed equally to this work and share first authorship