key: cord-0024528-dfyt8d28
authors: Wang, Enhao; Li, Yang; Ming, Ruijie; Wei, Jiahui; Du, Peiyu; Zhou, Peng; Zong, Shimin; Xiao, Hongjun
title: The Prognostic Value and Immune Landscapes of a m(6)A/m(5)C/m(1)A-Related LncRNAs Signature in Head and Neck Squamous Cell Carcinoma
date: 2021-11-30
journal: Front Cell Dev Biol
DOI: 10.3389/fcell.2021.718974
sha: bda357e7ad036ae5a71f209be14016a4e61aaf73
doc_id: 24528
cord_uid: dfyt8d28

Background: N6-methyladenosine (m(6)A), 5-methylcytosine (m(5)C) and N1-methyladenosine (m(1)A) are the main RNA methylation modifications involved in the progression of cancer. However, it is still unclear whether m(6)A/m(5)C/m(1)A-related long non-coding RNAs (lncRNAs) affect the prognosis of head and neck squamous cell carcinoma (HNSCC). Methods: We summarized 52 m(6)A/m(5)C/m(1)A-related genes, downloaded 44 normal samples and 501 HNSCC tumor samples with RNA-seq data and clinical information from The Cancer Genome Atlas (TCGA) database, and then searched for m(6)A/m(5)C/m(1)A-related genes co-expressed lncRNAs. We adopt the least absolute shrinkage and selection operator (LASSO) Cox regression to obtain m(6)A/m(5)C/m(1)A-related lncRNAs to construct a prognostic signature of HNSCC. Results: This prognostic signature is based on six m(6)A/m(5)C/m(1)A-related lncRNAs (AL035587.1, AC009121.3, AF131215.5, FMR1-IT1, AC106820.5, PTOV1-AS2). It was found that the high-risk subgroup has worse overall survival (OS) than the low-risk subgroup. Moreover, the results showed that most immune checkpoint genes were significantly different between the two risk groups (p < 0.05). Immunity microenvironment analysis showed that the contents of NK cell resting, macrophages M2, and neutrophils in samples of low-risk group were significantly lower than those of high-risk group (p < 0.05), while the contents of B cells navie, plasma cells, and T cells regulatory (Tregs) were on the contrary (p < 0.05). In addition, patients with high tumor mutational burden (TMB) had the worse overall survival than those with low tumor mutational burden. Conclusion: Our study elucidated how m(6)A/m(5)C/m(1)A-related lncRNAs are related to the prognosis, immune microenvironment, and TMB of HNSCC. In the future, these m(6)A/m(5)C/m(1)A-related lncRNAs may become a new choice for immunotherapy of HNSCC.

Head and neck squamous cell carcinoma (HNSCC) represents the sixth most common malignant tumor in the world . These tumors originate from the mucosa of various anatomical parts of the oral cavity, pharynx, larynx and sinonasal tract and are heterogeneous in etiology and clinicopathological characteristics, and these highly aggressive malignant tumors affect more than 830,000 patients worldwide each year Qin et al., 2021) . HNSCC is usually associated with exposure to tobacco-derived carcinogens, excessive alcohol consumption, or both. Tumors appearing in the oropharynx are increasingly related to human papillomavirus (HPV) carcinogenic strains, mainly HPV-16, and to a lesser extent HPV-18 and other strains of infection (Riva et al., 2021) . In current treatment methods, only about 50% of HNSCC patients survive for more than 5 years. Therefore, it is urgent to analyze the heterogeneity of HNSCC to choose different treatment options. In particular, there is an urgent need for reliable biomarkers with predictive and prognostic value (Rasheduzzaman et al., 2020) .

Some evidences show that epigenetic changes (including: DNA methylation, histone covalent modification, chromatin remodeling, non-coding RNA (ncRNA) activity and RNA chemical modifications) are often related to oral carcinogenesis, tumor malignancy and resistance to treatment (Castilho et al., 2017; Romanowska et al., 2020) . Increasing evidence suggests that dynamic RNA modifications pathways are also misregulated in human cancers (including HNSCC, lung cancer, breast cancer, pancreatic cancer, endometrial cancer, et al.) and may be ideal targets of cancer therapy (Barbieri and Kouzarides, 2020; Zhou et al., 2020; Huang et al., 2021) . Nowadays, over 100 types of chemical modifications have been identified in cellular RNAs, including messenger RNAs (mRNAs), transfer RNAs (tRNAs), ribosomal RNAs (rRNAs), and ncRNAs, most of these modifications contribute to pre-mRNA splicing, nuclear exporting, transcript stability and translation initiation of eukaryotic cellular RNAs. Among these modifications, methylation modification is the most abundant which include N6-methyladenosine (m 6 A), 5methylcytosine (m 5 C), N1-methyladenosine (m 1 A) (Shi et al., 2020) . M 6 A is the most abundant form of methylation modification in eukaryotic cellular RNAs, and it is also the most thoroughly studied type of RNA modification. The sequence near the m6A modification site on the mRNA is highly conservative, mainly occurring on the adenine of the RRACH (R: purine; A: m 6 A; H: non-guanine) sequence, and its function is determined by the methyltransferase (writers: methyltransferase like 3 (METTL3), methyltransferase like 14 (METTL14) and WT1 associated protein (WTAP), et al.), demethylase (erasers: fat mass and obesity associated protein (FTO) and alkB homologue 5 (ALKBH5), et al.) and binding protein (readers: YTH domain containing family, the insulin like growth factor 2 mRNA binding protein (IGF2BP) family and heterogeneous nuclear ribonucleoproteins A2/B1 (HNRNPA2B1), et al.) are jointly regulated . M 6 A modification is involved in a variety of biological processes, such as stem cell differentiation, cell division, gametogenesis and biological rhythms, et al., as well as the occurrence of a variety of diseases, including tumors, obesity and infertility Zhang M. et al., 2020; Song et al., 2020) . M 5 C modification is another prevalent RNA modification in multiple RNA species, including mRNAs, tRNAs, rRNAs, and ncRNAs . M 5 C modification is distributed with mRNA, enriched around 5′UTR and 3′UTR, and conserved in tRNA and rRNA. It is dynamically regulated by related enzymes, including methyltransferase (writers: NOL1/NOP2/SUN domain family member (NSUN) family and DNA methyltransferase like 2 (DNMT2), et al.), demethylases (erasers: tet methylcytosine dioxygenase 2 (TET2), et al.), and binding proteins (readers: Aly/ REF export factor (ALYREF), et al.) (Wiener and Schwartz, 2021) . Up to now, accumulated studies have shown that m 5 C is involved in various RNA metabolism, including mRNA export, RNA stability and translation. The absence of m 5 C modification in organisms may lead to mitochondrial dysfunction, stress response defects, gametogenesis and embryogenesis frustration, and abnormal development of nerves and brain, which is related to cell migration and tumorigenesis (Chellamuthu and Gray, 2020; Vissers et al., 2020) . The majority of these m 1 A sites are located within the 5′-untranslated region (5′ UTR), m 1 A modification has been identified in tRNA, rRNA, mRNA, long non-coding RNA (lncRNA) and mitochondrial (mt) transcripts, It can also be dynamically regulated by methyltransferases, demethylases and binding proteins (Zhang and Jia, 2018) .

LncRNAs are a class of RNA molecules with transcript length exceeding 200nt. it is generally believed that they do not encode proteins, but participate in the regulation of protein-encoding genes in many aspects (epigenetic regulation, transcriptional regulation and post-transcriptional regulation, etc.) in the form of RNA. More and more lncRNAs have been found to play a role in the initiation and development of HNSCC, suggesting that they can be used as novel biomarkers and therapeutic targets to provide more effective diagnosis, prognosis and treatment for HNSCC patients, such as HOTAIR, UCA1, FOXC1, AFAP1-AS1, et al. (Cossu et al., 2019; Wang S. et al., 2020; Wang et al., 2020c.) . Researchers have found that methylation modification of lncRNA regulates numerous processes affecting cancer cell activity by enhancing and increasing lncRNA stability, increasing lncRNA nuclear accumulation, and promoting decay of some lncRNAs (Dai et al., 2020) . However, there is currently no evidence to explain the correlation between three major RNA methylation modifications related lncRNAs and HNSCC. In this study, we integrated 44 normal samples and 501 tumor samples with clinical trials data from The Cancer Genome Atlas (TCGA). We used bioinformatics and statistical analysis methods to explore whether m 6 A/m 5 C/m 1 A-related lncRNAs were related to the diagnosis and prognosis of HNSCC.

RNA-seq data and clinical information of HNSCC patients were downloaded from The Cancer Genome Atlas (TCGA) database Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 2 (https://cancergenome.nih.gov/) (Wang et al., 2020d) , which contained data of 501 HNSCC and 44 non-tumor tissues. Raw count data was downloaded for further analysis. In addition, based on previous publications, the expression matrix of 52 m 6 A/ m 5 C/m 1 A-related genes was extracted from TCGA database, including expression data on writers, erasers and readers ( Table 1) . We obtain the gene annotation file of Genome Reference Consortium Human Build 38 (GRCh38) from the GENCODE website (https://www.gencodegenes.org/human/) for the annotation of the mRNAs in the TCGA database (Frankish et al., 2019) .

The simple nucleotide variation files with masked somatic mutation data were acquired from the Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov/) (Gong et al., 2018) .

We downloaded 33 types of tumor HTSeq-FPKM GDC Hub files, Phenotype GDC Hub files, survival data GDC Hub files from the UCSC Xena database (http://xena.ucsc.edu/) , and downloaded immune subtype files from the TCGA Pan-Cancer (PANCAN) database. We downloaded the nci 60 RNA seq files from Processed Data Set and Compound activity: DTP NCI-60 files from the CellMiner database (http://discover. nci.nih.gov/cellminer/home.do) (Reinhold et al., 2019) .

2.2.1. mRNAs-lncRNAs Co-Expression Network and Prognostic Value of m 6 A/m 5 C/m 1 A-Related lncRNAs After obtaining the expression matrix of 52 m 6 A/m 5 C/ m 1 A-related genes and expression matrix of lncRNAs from the TCGA database, we screened out the lncRNAs coexpressed with m 6 A/m 5 C/m 1 A-related genes by setting the absolute value of correlation coefficient > 0.5 and p value < 0.001. The univariate cox regression analysis was performed to screen out prognostic value of m 6 A/m 5 C/m 1 A-related lncRNAs.

Based on the expression level of these 44 lncRNAs with prognostic value in tumor samples, we used the R package "ConsensusClusterPlus" to classify tumor samples into different groups and to analyze the prognosis by cluster analysis of samples in the database . Principal Component Analysis (PCA) was applied to the typed samples to observe whether the samples could be intuitively distinguished (Rapa et al., 2021) . In order to understand the prognosis of samples with different classifications, the R package "survival" and "Survminer" (Zhang L.-P. et al., 2020) were used to analyze the prognosis of the samples after classifications. Clinical data were included to observe whether different clinical features differed in different classifications. The NMT classification comes from the TNM classification of malignant tumours (TNM) (https://www.uicc.org/resources/tnm), which is a globally recognised standard for classifying the extent of spread of cancer. Each individual aspect of TNM is termed as a category: T category describes the primary tumor site and size; N category describes the regional lymph node involvement; M category describes the presence or otherwise of distant metastatic spread. Gene enrichment analysis (GSEA) was used to understand the differentially functional enrichment pathways of samples with different types (Tang et al., 2021) . In order to explore the characteristics of immune microenvironment of different types of samples, the R package "limma" was used to filter the differences in immune cell infiltration of different types of samples . In addition, we further analyzed the expression of 29 immune checkpoint moleclues in different types of samples, including CD70, TNFSF9, FGL1, CD276, NT5E, HHLA2, VTCN1, TNFRSF18, CD274, PDCDL1LG2, IDO1, CTLA4, ICOS, HAVCR, PDCD1, LAG3, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, VSIR, ICOSLG, TNFSF14, TMIGD2, CD27, NCR3, BTLA, CD40LG, CD40, TNFRSF4 (Zhang C. et al., 2020) . Moreover, we also analyzed the differences in the expression levels of the tumor suppressor gene TP53 among the subtypes.

The least absolute shrinkage and selection operator (LASSO) cox regression was used to select the most significant prognostic lncRNAs in 44 m6A/m5C/m1A-related lncRNAs for predictive model construction (Wang S. et al., 2020) . Lasso is a linear regression method based on L1-regularization. L1regularization can reduce the complexity of the model and the risk of over-fitting of the model. The explanation in regression analysis is that if two characteristic variables are highly correlated, the regression analysis without regularization may give them approximately equal weight, so the weighted result is still small, but due to the large weight, the over-fitting problem is caused. L1regularization greatly increases the stability of modeling and achieves the purpose of selecting characteristic variables. The most outstanding advantage of Lasso lies in the penalized regression analysis of all variable coefficients, the relatively unimportant independent variable coefficients become 0, which is excluded from modeling (Vaid et al., 2021; Yu et al., 2021; Hoq et al., 2021) . However, the 44 lncRNAs are too large for constructing a prognosis model, so we used LASSO regression analysis to further compress the 44 lncRNAs to get the lncRNAs with the most prognostic value, and their correlation coefficients.

Combining respective correlation coefficients of obtained m 6 A/ m 5 C/m 1 A-related lncRNAs, the respective gene expression level of these lncRNAs in each tumor sample, and the risk score calculating formula (The risk score calculating formula is:

where Coefi means the coefficients, xi is the FPKM value of each m 6 A/m 5 C/ m 1 A-related lncRNAs), we can grade each sample. All samples can be sorted according to their respective scores. At last, with the median of sample scores as the boundary, tumor samples with scores lower than the median of all scores were divided into low risk group, while tumor samples with scores higher than the median of all scores were divided into high risk group. Finally, a prognostic model of risk score was constructed based on these obtained m 6 A/m 5 C/m 1 A-related lncRNAs. The R package "survival" was used to analyze the survival of the high-risk group and low-risk group in HNSCC. Tumor samples were then divided into a training set (249 samples) and a testing set (249 samples), and there was no statistically significant difference in clinical characteristics between the two groups (Supplementary Table S1 ). The survival conditions of the high-risk group and lowrisk group were analyzed in the training set, and the results were verified in the testing set. The univariate and multivariate regressions were further used to analyze whether risk score in the training set and validation set could be used as an independent prognostic indicator. In addition, PCA and T-Distributed Neighbor Embedding (T-SNE) analysis were used to visualize whether the high-risk group and low-risk group could be distinguished by risk score. Receiver operating characteristic (ROC) curve was used to test model effectiveness (Longlong et al., 2020) . Next, we evaluated the clinical characteristics of the high-risk group and the lowrisk group: age, gender, HPV-status, tumor stage, tumor grade, and tumor TNM, and performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of the genes that differed between the two groups .

Considering the strength expression of markers in samples. If it is a weak expression, It might be hard to be detected in clinical diagnosis. Therefore, we searched these filtered lncRNAs included in this study in TIMER: Tumor IMmune Estimation Resource (http://timer.cistrome.org/) .

In order to explore the immune landscapes of the high-risk group and low-risk group, we analyzed the correlation between risk score and immune cells by Spearman analysis. The R Package "Limma" was used to filter the immune cell infiltration of samples with different risk scores. In addition, we analyzed the expression levels of 29 immune checkpoint moleclues in samples from different risk groups, including CD70, TNFSF9, FGL1, CD276,   NT5E, HHLA2, VTCN1, TNFRSF18, CD274, PDCDL1LG2,  IDO1, CTLA4, ICOS, HAVCR, PDCD1, LAG3, SIGLEC15,  TNFSF4, TNFRSF9, ENTPD1, VSIR, ICOSLG, TNFSF14,  TMIGD2, CD27, NCR3 , BTLA, CD40LG, CD40, TNFRSF4.

Based on the tumor mutational burden (TMB) information of each sample (Bonfield et al., 2021) , combined with the survival information analysis of patients, all samples were divided into the high-TMB group and the low-TMB group, and the survival of the two groups was analyzed. The survival of the high-TMB + high-risk group, the high-TMB + low-risk group, the low-TMB + high-risk group and the low-TMB + low-risk score group were jointly analyzed to observe the influence of the risk score and TMB on the survival of the patients. Genes with significant mutation rates in the high-risk group and low-risk group were further analyzed to find the ones with the most significant mutation rates for further analysis.

Based on the expression levels of each independent lncRNA in the 6 m 6 A/m 5 C/m 1 A-related lncRNAs Prognostic signatures, the patient samples were divided into high expression group and low expression group. The R package "survival" was used to obtain survival curves for each independent lncRNA. There are six types of immune infiltration in human cancers, including Wound Healing (Immune C1), IFN-gamma Dominant (Immune C2), Inflammatory (Immune C3), Lymphocyte Depleted (Immune C4), Immunologically Quiet (Immune C5) and TGF-beta Dominant (Immune C6) (Soldevilla et al., 2019; Zhao et al., 2021) . Based on the relevant files downloaded from the UCSC Xena database, the R package "CIBERSORT" and "ESTIMATE" (Wang S. et al., 2020) were used to obtain relevant data (Stromal Score, Immune Score, ESTIMAT Score) of tumor purity. Furthermore, the correlation between the expression levels of 6 m 6 A/m 5 C/m 1 A-related lncRNAs and tumor purity in the immune microenvironment of HNSCC was further analyzed.

The R package "Prophetic" was used to analyze the drug sensitivity analysis of samples from different risk groups (Ding et al., 2021) . Combined with the relevant data downloaded from the CellMiner Database, the correlation analysis of drug sensitivity and the expression level of 6 m 6 A/m 5 C/m 1 A-related lncRNAs was further carried out .

The flow chart of our work was shown in Figure 1 .

Using the downloaded file from the "GENCODE" website, we obtained 14,086 lncRNAs expression matrix. In addition, 52 expression matrix of m 6 A/m 5 C/m 1 A-related genes were screened out from the expression matrix of mRNAs. Through co-expression analysis, we screened out 147 lncRNAs co-Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 expressed with m 6 A/m 5 C/m 1 A-related genes ( Figure 2A) . We further screened out 44 lncRNAs with potential prognostic value from 147 co-expressed lncRNAs. The forest plot showed that all 44 lncRNAs were protective factors (hazard ratio < 1) in patients with HNSCC ( Figure 2B) . 

By cluster analysis, the tumor samples in the TCGA database were classified into two types, consensus matrix for optimal K 2 ( Figure 3A) . PCA was used to represent the characteristics of the typing samples, and it was found that the samples could be visually distinguished ( Figure 3B ). The Kaplan-Meier survival curves showed that the prognosis of patients in Cluster 2 was significantly better than that of patients in Cluster 1 (p 0.03) ( Figure 3C ), and the two clinical characteristics of patients' HPV status and tumor grade were significantly different among different types (p < 0.001) ( Figure 3D ). GSEA analysis showed that the pathways of Cluster 2 were mainly enriched in CELL_CYCLE, LYSINE_DEGRADATION, P53_SIGNALING_PATHWAY, et al.

The functional enrichment pathways of Cluster1 were mainly enriched in ECM_receptor_interaction, FOCAL_ADHESION, OCAL_ADHESION, et al. ( Figure 3E ). According to the immune microenvironments of the two subtypes, we can found the infiltration of T cells follicular helper and T cells regulatory in Cluster 1 significantly less than Cluster 2 (p 0.005; p 0.009), while Macrophages M0 infiltrate significantly more than Cluster 2 (p 0.05) ( Figure 3F ). We found that the expression levels of 9 immune checkpoint genes were significantly different between the two subtypes (p < 0.05). The expression levels of NT5E, ICOS, PDCD1LG2, HAVCR2, TNFSF4, TNFRSF9, EntPD1 and TNFSF14 were higher in Cluster 1, and TNFRSF18 was higher in Cluster 2 ( Figure 3G ). In addition, the expression level of tumor suppressor gene TP53 in Cluster 2 was significantly higher than that in Cluster 1 (p < 0.001) ( Figure 3H ).

We got 44 lncRNAs with prognostic value, but it is too difficult to use all lncRNAs to build a prognostic model, so we compressed 44 lncRNAs with Lasso analysis method, and finally got 6 lncRNAs with the most prognostic value to build our prognostic model: AL035587.1, AC009121.3, AF131215.5, FMR1-IT1, AC106820.5, PTOV1-AS2 ( Figures 4A,B) . The correlation coefficient was shown in Table 2 .

The survival time and survival rate of patients in the high-risk subgroup were significantly lower than those in the low-risk subgroup (p 5.262e-04) ( Figure 4C) , which was also verified in samples from the training set and the testing set (Train group: p 0.008; Test group: p 0.02) (Figures 4D,E) . It was found in the training set and testing set that the prognosis of patients in the high-risk subgroup was worse than that in the low-risk subgroup, and the expression levels of 6 lncRNAs in the lowrisk subgroup were significantly higher than that in the highrisk subgroup (Figures 4F,G) . The univariate and multivariate regression analysis results showed that risk score, as an independent prognostic indicator, had a certain predictive value relative to individual and all clinical characteristics (age, gender, tumor grade, tumor stage) in both the training set and the testing set (p < 0.01) (Supplementary Figures  S1A-D) . The results of PCA and T-SNE show that the use of risk score can intuitively distinguish the high-risk group from the low-risk group in the training set and the testing set (Supplementary Figures  S1E-H) . As shown in Supplementary Figure S2 , the area under curve (AUC) at 5 years of overall survival for the LncSig (our study) is 0.607, which was significantly higher than that of JianglncSig (AUC 0.602) and WanglncSig (AUC 0.589). These results indicate that LncSig has better prognostic performance in predicting survival than the two recently published lncRNA signatures (Wang et al., 2018; Jiang et al., 2020) .

In combination with clinical features, analysis showed that patients in the tumor grade G3-4, N0, N1-3, stage I-II, stage Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 5 III-IV, T1-2, T3-4, the survival of high-risk subgroup was significantly worse than low-risk subgroup (p < 0.05) ( Figures  5A-H) . The heatmap showed the expression of 6 selected lncRNAs in patients with high and low risk. We observed the statistically significant differences between the high-risk and lowrisk groups in gender, HPV-status, tumor stage and tumor grade (p < 0.05) ( Figure 5I ). The results of sankey diagram showed that most of the patients in Cluster 2 belonged to the low-risk group and had a higher survival rate than Cluster 1, which also supported the results of survival status analysis of the above classifications ( Figure 5J ). GO analysis showed that cell functions enriched by different genes between the high-risk group and the low-risk group were as follows: B cell mediated immunity, complement activation, complement activation classical pathway, humoral immune response mediated by circulating immunoglobulin, immunoglobulin mediated immune response ( Figure 5K) . KEGG analysis showed that the signal pathways enriched by differential genes between the high-risk group and the low-risk group were: drug metabolism-cytochrome P450, estrogen signaling pathway, leukocyte transendothelial migration, retinol metabolism, signaling pathway regulating pluripotency of stem cells ( Figure 5L) .

The gene expression box plots (Supplementary Figure S3 ) of FMR1-IT1 and PTOV1-AS2 were found in the Module of exploring associations between gene expression and tumor features in TCGA from TIMER, which showed that the expression of FMR1-IT1 and PTOV1-AS2 in HSNCC tissues is significantly higher than that in normal tissues.

We found that 12 immune checkpoint genes were significantly different between the two risk groups (p < 0.05). CD70, ICOS, HAVCR2, PDCD1, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 6 TNFSF14, CD27, BTLA, CD40LG. The expression of BTLA was higher in low-risk group. CD70, ICOS, HAVCR2, PDCD1, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, TNFSF14, CD27 and CD40LG were higher expressed in high-risk group ( Figure 6A ). Among them, the overlapping immune checkpoint genes ICOS, HAVCR2, TNFSF4, TNFRSF9, ENTPD1 and TNFSF14 were all increased in Cluster 1 and high-risk subgroups.

Immunity microenvironment analysis showed that the contents of NK cell resting, macrophages M2, and neutrophils in samples of low-risk group were significantly lower than those of high-risk group (p < 0.05), while the contents of B cells navie, plasma cells, and T cells regulatory (Tregs) were on the contrary (p < 0.05) ( Figure 6B ). The correlation analysis between risk score and immune cell infiltration showed that neutrophils, dendritic cells activated, NK cells resting, dendritic cells resting, macrophages M2, mast cells activated were positively correlated with risk score. B cells naive, plasma cells, T cells helper and T cells regulatory (Tregs) were negatively correlated with risk score (p < 0.05) (Figures 6C-L) .

TMB analysis showed that the survival of the high-TMB group was significantly lower than that of the low-TMB group (p 0.002) Figure 7A ). Combined analysis of risk score showed that patients in the low-TMB + low-risk score group had a significantly better prognosis than those in the other three groups (p < 0.001), suggesting that patients with low TMB tend to have a better prognosis ( Figure 7B ). Looking for the significant mutant genes in the high-risk group and the low-risk group, it was found that the mutation rate of TP53 was the highest in both groups, and the expression of TP53 in the high-risk group was significantly lower than that in the low-risk group (p 1.3e-07) ( Figures 7C-E) . The combined analysis of TP53 and risk score showed that the expression level of TP53 in the high-risk group was higher than that in the low-risk group (p 1.3e-07). Patients in the low-risk subgroup of TP53 wild-type had significantly better outcomes than those in the other three groups (p 0.002) ( Figure 7F ).

The Kaplan-Meier survival curves showed that there was no significant difference between high-expression and lowexpression of AF131215.5. However, the high-expression of AL035587.1, AC009121.3, FMR1-IT1, AC106820.5, PTOV1-AS2 had the better overall survival than low-expression in TCGA database (p < 0.05) ( Figures 8A-F) .

According to the results of the immunoinfiltration analysis, the high expression of AL035587.1 was considered to be related to the type of C3 and C6 infiltration according to the results of Figure 9A . The high expression of AC106820.5 was considered to be correlated with the type of C3 and C6 infiltration (p < 0.01). Moreover, the tumor microenvironment (TME) was assessed by estimation analysis. The expression levels of PTOV1-AS2, AC009121.3, FMR1-IT1 and AC106820.5 were negatively correlated with the stromalscore ( Figure 9B ). In addition, AL035587.1 and AF131215.5 were positively correlated with immunescore, while AC106820.5 was negatively correlated. Estimatescore (combining stromalscore and immunescore) showed that the expression levels of AC009121.3, AC106820.5 and FMR1-IT1 were negatively correlated with estimatescore, while AL035587.1 and AF131215.5 were positively correlated ( Figure 9B ) (p < 0.05).

Drug sensitivity analysis showed that samples from the low-risk group had higher sensitivity to Paclitaxel, Docetaxel and Gefitinib and lower sensitivity to Methotrexate than those from the highrisk group (Figures 10A-D) (p < 0.001). Two m 6 A/m 5 C/ m 1 A-related lncRNAs: PTOV1-AS2 and FMR1-IT1 were extracted from the CellMiner Database NCI 60 RNA seq and Compound activity: DTP NCI-60. It was found that the sensitivity of drug Nelarabine had the highest correlation with the expression level of PTOV1-AS2 (Correlation 0.418, p < 0.001), and the sensitivity of Nelarabine was also the highest associated with FMR1-IT1 expression (Correlation 0.337, p 0.008) ( Figure 10E ). 

Studies have shown that dynamic RNA methylation, the modification events (such as m 6 A, m 5 C, and m 1 A) are involved in the progression of tumors, including promoting cancer cell proliferation or regulating invasiveness and metastatic potential (Barbieri and Kouzarides, 2020) . m 6 A can also be involved in tumor proliferation, migration, drug resistance, or as a prognostic marker of tumor Wang M. et al., 2020; Paramasivam and Priyadharsini, 2021; Xu et al., 2021) . Some studies have shown that there is an interaction between RNA methylation and lncRNA in tumors, and m 6 A reader YTHDF2mediated degradation of lncRNA FENDRR promoted cell proliferation in endometrioid endometrial carcinoma (Shen Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 9 et al., 2021). m 6 A-mediated upregulation of LINC00857 promoted pancreatic cancer tumorigenesis (Meng L. et al., 2021) . LncRNA KB-1980E6.3 maintained breast cancer stem cell stemness via interacting with IGF2BP1 (Zhu P. et al., 2021) . Abnormal m 5 C modification of lncRNA H19 mediated by NSUN2 was associated with poor differentiation of hepatocellular carcinoma (Sun et al., 2020) . Conversely, lncRNAs can also affect the function of genes related to RNA methylation modifications. LncRNA MALAT1 contributed to thyroid cancer progression through the upregulation of IGF2BP2 by binding to miR-204 (Ye et al., 2021) . LncRNA LINRIS can stabilize IGF2BP2 and promote aerobic glycolysis in colorectal cancer . The above evidences indicate the interaction between RNA methylation modifications and lncRNAs in the tumorigenesis of a variety of cancers. However, the evidence of how the RNA methylation modification involved in the development of HNSCC by influencing the lncRNA is still insufficient. Therefore, this study aims to explore whether the relevant lncRNAs are involved in the development of the HNSCC.

We downloaded the data of 44 Normal samples and 501 tumor samples from the TCGA database to search for m 6 A/m 5 C/ m 1 A-related lncRNAs and we found 44 lncRNAs with prognostic value. Tumor samples were divided into Cluster 1 and Cluster 2 according to the expression levels of these 44 lncRNAs in all tumor samples, and the survival of Cluster 2 was significantly better than that of Cluster 1. We further used the lasso cox regression to construct the prognostic signature using 6 m 6 A/m 5 C/m 1 A-related lncRNAs to predict survival in patients with HNSCC. Based on the correlation coefficients of the 6 m 6 A/ m 5 C/m 1 A-related lncRNAs and their expression levels in each tumor sample, we obtained risk scores for each tumor sample and constructed a prognostic signature. According to Figure 5J the sankey diagram of the relationship among the cluster, risk score and survival state, it can be found that most patients in high risk group come from cluster 1, while most patients in cluster 2 are assigned to low risk group. The survival time of patients in the high risk group was significantly shorter than that in the low risk group (p < 0.05), which also proved that both the clustering type Top 20 genes with the most significant mutation rates in high-risk subgroup and low-risk subgroup. (E) The differential expression level TP53 between high-risk group and low-risk group. (F) Kaplan-Meier curve analysis of overall survival was shown for patients classified according to TP53 mutation status and risk score. TP53-wt: wild type of TP53 sequence; TP53-mt: mutation type of TP53 sequence.

Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 11 composed of 44 lncRNAs and the prognosis model composed of 6 lncRNAs can well predict the prognosis of tumor patients. Differences in survival were also demonstrated in the training and validation sets of randomly assigned tumor samples. The univariate and multivariate cox regression analysis showed that risk score could be used as a prognostic marker for HNSCC independent of age, sex, tumor grade, and tumor stage. This suggested that the prognostic signature constructed by 6 lncRNAs had certain prognostic value. These results supply a reference for further functions and mechanisms of prognostic signatures in HNSCC.

Six m 6 A/m 5 C/m 1 A-related lncRNAs were obtained, which were AL035587.1, AC009121.3, AF131215.5, FMR1-IT1, AC106820.5 and PTOV1-AS2. The results in Supplementary Figure S3 showed that FMR1-IT1 and PTOV1-AS2 can be detected in clinical diagnosis. AF131215.5 is one of lncRNAs in lung adenocarcinoma (LUAD)-related competing endogenous RNA networks, and showed an independent prognostic value of overall survival for patients with LUAD (Hou and Yao, 2021) . FMR1 intronic transcript 1 (FMR1-IT1) is also known as FMR1 intronic transcript 1 (non-protein coding) or Fragile X Mental Retardation Syndrome Protein (FMRP) or FMR1, there were evidence that FMRP regulated tumor invasiveness-related pathways in melanoma cells by impacting cell migration, invasion and adhesion (Zalfa et al., 2017) , and its overexpression was associated with lung metastasis of murine breast cancer (Lucá et al., 2013) . However, there is still not enough evidence to clarify the relationship between tumorigenesis and prognosis of these 6 lncRNAs and HNSCC. We hope that our work can explain the role of these 6 lncRNAs in the carcinogenesis and development of HNSCC and clarify their prognostic value.

Based on the expression levels of these 6 m 6 A/m 5 C/ m 1 A-related lncRNAs in the samples, we further studied the tumor microenvironment of the samples and analyzed the expression levels of 29 immune checkpoint genes. There were significant differences in the expression of various immune cells and immune checkpoint genes between the high-risk subgroup and the low-risk subgroup. Among them, the overlapping immune checkpoint genes of Cluster 1 and high-risk subgroups were ICOS, HAVCR2, TNFSF4, TNFRSF9, ENTPD1 and TNFSF14. Their expression levels were all elevated, revealing the characteristics of tumor immune microenvironment between different subgroups. Currently, only a small number of HNSCC patients benefit from immunotherapy, and the need to discover new biomarkers to optimize treatment strategies has become increasingly important in clinical practice. With the rise of the application of immune Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 checkpoint inhibitor in tumor treatment, such as Programmed Cell death-1 (PD-1) inhibitors (Kitamura et al., 2020) , there are many ongoing trials in HNSCC that focus on identifying new biomarkers (Gavrielatou et al., 2020) . As a member of the B7 superfamily, ICOS is a co-stimulatory receptor for T-cell enhancement. In preclinical studies, ICOS agonist monoclonal antibodies have been shown to enhance inhibitory checkpoint blockade. On the contrary, the antagonistic monoclonal antibody against ICOS can not only inhibit the lymphoid tumor cells expressing ICOS, but also inhibit the immunosuppressive Treg (Amatore et al., 2018) . ENTPD1 (CD39) converts extracellular ATP into adenosine, thereby inhibiting T cell effector functions through the adenosine receptor A2A. The CD39 inhibitor POM-1 inhibits adenosine production and reduces T cell suppression. The adenosine pathway can be targeted to multiple myeloma, and blocking this pathway can replace PD1/PDL1 to inhibit multiple myeloma and other blood system cancers . At present, TNFSF14-based therapy has shown great efficacy in cancer immunotherapy, which can modify tumor microenvironment by normalizing tumor vessels and significantly improve the infiltration of effector infiltrating lymphocytes, thus reducing tumor burden and forming lasting FIGURE 9 | Association of expression of prognostic lncRNAs with with immune response, TME. *, **, *** represent p < 0.05, p < 0.01 and p < 0.0001, respectively.

(A) Association of prognostic lncRNAs expression with immune subtypes in HNSCC patients. (B) Correlation of expression of prognostic lncRNAs with TME (Stromal score, Immune score, and ESTIMATE score).

Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 718974 13 anti-tumor memory (Skeate et al., 2020) . These evidences also suggest the possibility of selecting differential immune checkpoint genes as therapeutic targets for HNSCC. In addition, we have analyzed the joint tumor mutation burden and risk score model. It was found that the prognosis of patients with high burden of tumor mutations was worse than that of patients with low burden, while the prognosis of patients with high burden of mutations in the high-risk subgroup was the worst, and the prognosis of patients with low burden of mutations in the low-risk subgroup was the best. This undoubtedly indicated the importance of tumor mutation burden in predicting prognosis. Moreover, we observed that the tumor suppressor gene TP53 was highly expressed in Cluster 2 and the low-risk subgroup. It is generally believed that TP53 gene can be divided into wild type and mutant type. The wild-type TP53 gene, namely tumor suppressor gene, has a control and negative regulation effect on cell proliferation. However, the mutant TP53 gene lost its tumor suppressor effect and promoted the transformation of cells to malignancy (Duffy et al., 2020; Zhu et al., 2020; Meng X. et al., 2021) . Our results showed that TP53 mutant patients in the high-risk subgroup had the worst prognosis, and TP53 wild-type patients in the low-risk subgroup had the best prognosis, which also indicated that higher mutation rate of TP53 was associated with poorer prognosis. In addition, some studies have revealed the connection between lncRNA and TP53. Huarte, M, et al. found that LincRNA-p21 is the first lncRNAs identified as being induced by wt-p53 transcription. It acts as a transcriptional repressor in the p53 pathway (Huarte et al., 2010) . Chang et al. (2011) found that lncRNA XIST and the p53-inducible lncRNA NEAT1 act as ceRNA for miR-34a, which leads to the activation of epithelial-mesenchymal transition. The research of Olivero et al. (2020) showed that TP53 activates the lncRNA pvt1b to inhibit myelocytomasis and suppress tumorigenesis. It is worth noting that targeting mut-TP53/ncRNA signals may have a great potential impact on the treatment of cancer patients. Interestingly, lncRNA and immune checkpoint are also closely related, lncRNAs are involved in the regulation of various cellular pathways and immune checkpoints (Di Martino et al., 2021) . Tang et al. (2017) found that lncRNA AFAP1-AS1 was coexpressed with PD-1 lymphocytes associated with nasopharyngeal carcinoma, which resulted in invasive and poor prognosis phenotype. In hepatocellular carcinoma, lncRNA Tim3 mediates T cell exhaustion and inhibits T cell immune response, which leads to tumor immunosuppression through TIM-3 specific binding (Ji et al., 2018) . In the future, the potential combination of lncRNAs targeted drugs and tumor immune checkpoint inhibitor antibodies may be used for experimental treatment of human cancer. Considering that the six m 6 A/m 5 C/ m 1 A-related lncRNAs were all protective factors of HNSCC (HR < 1), we analyzed the survival of the high-risk and low-risk subgroups of these 6 lncRNAs, and found that the prognosis of the high-expression subgroup of 5 lncRNAs was significantly better than that of the low-expression subgroup, which also proved their role as protective factors. The high expression of AL035587.1 and AC106820.5 is thought to correlate with the type of Immune C3 and C6 invasion, indicating their tumor suppressive effect. Estimatedscores (a combination of stromalscore and immunescore) showed that high expression levels of AC009121.3, AC106820.5 and FMR1-IT1 were significantly associated with lower tumor purity, while AL035587.1 and AF131215.5 had the opposite effect (Sawatani et al., 2020) .

Drug sensitivity analysis showed that Paclitaxel, Docetaxel, and Gefitinib had higher sensitivity in the low-risk group than in the high-risk group, while Methotrexate had higher sensitivity in the high-risk group. These four drugs are commonly used in the clinical treatment of HNSCC (Yin et al., 2019; Ham et al., 2020; Guigay et al., 2021) , and our results also demonstrate their value in the application of HNSCC. For the therapeutic potential of these 6 lncRNAs, we analyzed their sensitivity to different small molecule drugs. The results showed that PTOV1-AS2 was sensitive to four small molecular agents (Nelarabine, Decitabine, Acrichine, Dasatinib) (p < 0.05), and FMR1-IT1 was sensitive to five small molecular agents (Nelarabine, Chelerythrine, Tamoxifen, Tanespimycin, Sunitinib (p < 0.05). These eight drugs have certain anticancer activity and can be considered as candidates for the targeted therapy of HNSCC. Nelarabine is a purine analogue approved for the treatment of patients with T-cell lymphoblastic lymphoma (Teachey and O'Connor, 2020) . Decitabine is a DNA demethylation agent, which plays an unmethylated role in inducing apoptosis and activating autophagy in myeloid leukemia . Dasatinib is an antagonist of the SRC family of kinases, and may potentially induce T-cell response leading to improved survival of malignant pleural mesothelioma . Chelerythrine, Tamoxifen, Tanespimycin and Sunitinib have also been shown to have anticancer effects in a variety of cancers (Heng and Cheah, 2020; Wang and Liu, 2020; Ocadlikova et al., 2021; Su et al., 2021) . These evidences also indicate the prospect of targeting these two lncRNAs for the treatment of HNSCC.

The poor prognosis of the HNSCC affects the health of tens of millions of people every year, and improving the prognosis of the HNSCC is particularly important. Studies have shown that m 6 A/ m 5 C/m 1 A RNA methylation modifications was involved in the progression of cancer. Currently, many studies have chosen to use lncRNAs to establish prognostic markers, in order to provide new targets for diagnosis and treatment of cancers. Our study elucidates how m 6 A/m 5 C/m 1 A -related lncRNAs are related to the prognosis, immune microenvironment, and tumor mutation load of HNSCC. The prognostic signature can be used as an independent factor to better predict the prognosis of HNSCC, which may potentially become a new option for immunotherapy of HNSCC in the future.

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. 

Inducible Co-stimulator (ICOS) as a Potential Therapeutic Target for Anti-cancer Therapy

Role of RNA Modifications in Cancer

HTSlib: C Library for reading/writing High-Throughput Sequencing Data

Epigenetic Modifications and Head and Neck Cancer: Implications for Tumor Progression and Resistance to Therapy

p53 Regulates Epithelial-Mesenchymal Transition and Stem Cell Properties through Modulating miRNAs

The RNA Methyltransferase NSUN2 and its Potential Roles in Cancer

Evolution of Genomic and T-Cell Repertoire Heterogeneity of Malignant Pleural Mesothelioma under Dasatinib Treatment

Dynamic Transcriptomic M 5 C and its Regulatory Role in RNA Processing

Long Non-coding RNAs as Important Biomarkers in Laryngeal Cancer and Other Head and Neck Tumours

Crosstalk between RNA m6A Modification and Non-coding RNA Contributes to Cancer Growth and Progression

miRNAs and lncRNAs as Novel Therapeutic Targets to Improve Cancer Immunotherapy

Characterization of the Fatty Acid Metabolism in Colorectal Cancer to Guide Clinical Therapy

Targeting P53 for the Treatment of Cancer

GENCODE Reference Annotation for the Human and Mouse Genomes

Biomarkers for Immunotherapy Response in Head and Neck Cancer

Tumor Microenvironment Properties Are Associated with Low CD68-Positive Cell Infiltration and Favorable Disease-free Survival in EGFR-Mutant Lung Adenocarcinoma

Cetuximab, Docetaxel, and Cisplatin versus Platinum, Fluorouracil, and Cetuximab as First-Line Treatment in Patients with Recurrent or Metastatic Head and Neck Squamous-Cell Carcinoma (GORTEC 2014-01 TPExtreme): a Multicentre, Open-Label, Randomised, Phase 2 Trial

Methotrexate Plus or Minus Cetuximab as First-line Treatment in a Recurrent or Metastatic (R/M) Squamous Cell Carcinoma Population of the Head and Neck (SCCHN), Unfit for Cisplatin Combination Treatment, a Phase Ib-randomized Phase II Study Commence

Promote Tumor Proliferation of Bladder Cancer by Accelerating Pri-miR221/ 222 Maturation in m6A-dependent Manner

Chelerythrine Chloride Downregulates β-Catenin and Inhibits Stem Cell Properties of Non-small Cell Lung Carcinoma

Vocal Feature Extraction-Based Artificial Intelligent Model for Parkinson's Disease Detection

Potential Prognostic Biomarkers of Lung Adenocarcinoma Based on Bioinformatic Analysis

The Role of RNA N6-Methyladenosine Methyltransferase in Cancers

A Large Intergenic Noncoding RNA Induced by P53 Mediates Global Gene Repression in the P53 Response

Long Non-coding RNA Lnc-Tim3 Exacerbates CD8 T Cell Exhaustion via Binding to Tim-3 and Inducing Nuclear Translocation of Bat3 in HCC

A Novel Three-lncRNA Signature Predicts the Overall Survival of HNSCC Patients

Head and Neck Squamous Cell Carcinoma

Current Trends and Future Prospects of Molecular Targeted Therapy in Head and Neck Squamous Cell Carcinoma

Decitabine Downregulates TIGAR to Induce Apoptosis and Autophagy in Myeloid Leukemia Cells

The Role of N6-Methyladenosine mRNA in the Tumor Microenvironment

Identification and Analysis of Potential Key Genes Associated with Hepatocellular Carcinoma Based on

Predictive Value of the Texture Analysis of Enhanced Computed Tomographic Images for Preoperative Pancreatic Carcinoma Differentiation

The Fragile X Protein Binds M RNA S Involved in Cancer Progression and Modulates Metastasis Formation

CCR4, CCR8, and P2RY14 as Prognostic Factors in Head and Neck Squamous Cell Carcinoma Are Involved in the Remodeling of the Tumor Microenvironment

m6A-Mediated Upregulation of LINC00857 Promotes Pancreatic Cancer Tumorigenesis by Regulating the miR-150-5p/E2F3

Sunitinib Exerts In Vitro Immunomodulatory Activity on Sarcomas via Dendritic Cells and Synergizes with PD-1 Blockade

p53 Activates the Long Noncoding RNA Pvt1b to Inhibit Myc and Suppress Tumorigenesis

RNA N6-Methyladenosine: a New Player in Autophagy-Mediated Anti-cancer Drug Resistance

Tumor Microenvironment and Immune-Related Therapies of Head and Neck Squamous Cell Carcinoma

Bioactive Compounds in Cherry Tomatoes (Solanum Lycopersicum Var. Cerasiforme): Cultivation Techniques Classification by Multivariate Analysis

Protein Glycosylation in Head and Neck Cancers: From Diagnosis to Treatment

RNA Sequencing of the NCI-60: Integration into CellMiner and CellMiner CDB

HPV Meets APOBEC: New Players in Head and Neck Cancer

Head and Neck Squamous Cell Carcinoma: Epigenetic Landscape

Paclitaxel Potentiates the Anticancer Effect of Cetuximab by Enhancing Antibody-dependent Cellular Cytotoxicity on Oral Squamous Cell Carcinoma Cells In Vitro

Reader YTHDF2-Mediated Long Noncoding RNA FENDRR Degradation Promotes Cell Proliferation in Endometrioid Endometrial Carcinoma

Novel Insight into the Regulatory Roles of Diverse RNA Modifications: Re-defining the Bridge between Transcription and Translation

TNFSF14: LIGHTing the Way for Effective Cancer Immunotherapy

The Correlation between Immune Subtypes and Consensus Molecular Subtypes in Colorectal Cancer Identifies Novel Tumour Microenvironment Profiles, with Prognostic and Therapeutic Implications

Novel Insights into Adipogenesis from the Perspective of Transcriptional and RNA N6-Methyladenosine-Mediated Post-Transcriptional Regulation

High-dose Tamoxifen in High-Hormone-Receptor-Expressing Advanced Breast Cancer Patients: a Phase II Pilot Study

Aberrant NSUN2-Mediated m5C Modification of H19 lncRNA Is Associated with Poor Differentiation of Hepatocellular Carcinoma

Upregulation of SKA3 Enhances Cell Proliferation and Correlates with Poor Prognosis in Hepatocellular Carcinoma

Co-expression of AFAP1-AS1 and PD-1 Predicts Poor Prognosis in Nasopharyngeal Carcinoma

How I Treat Newly Diagnosed T-Cell Acute Lymphoblastic Leukemia and T-Cell Lymphoblastic Lymphoma in Children

Federated Learning of Electronic Health Records to Improve Mortality Prediction in Hospitalized Patients with COVID-19: Machine Learning Approach

The Epitranscriptome in Stem Cell Biology and Neural Development

m6A RNA Methylation Regulators Act as Potential Prognostic Biomarkers in Lung Adenocarcinoma

Upregulation of METTL14 Mediates the Elevation of PERP mRNA N6 Adenosine Methylation Promoting the Growth and Metastasis of Pancreatic Cancer

A Three-lncRNA Expression Signature Predicts Survival in Head and Neck Squamous Cell Carcinoma (HNSCC)

VDAC Upregulation and αTAT1-mediated α-tubulin A-cetylation C-ontribute to T-anespimycin-induced A-poptosis in Calu-1 C-ells

An Eight-CircRNA Assessment Model for Predicting Biochemical Recurrence in Prostate Cancer

RNA Sequencing Reveals the Expression Profiles of circRNA and Identifies a Four-circRNA Signature Acts as a Prognostic Marker in Esophageal Squamous Cell Carcinoma

GBX2, as a Tumor Promoter in Lung Adenocarcinoma, Enhances Cells Viability, Invasion and Migration by Regulating the AKT/ERK Signaling Pathway

LncRNA LINRIS Stabilizes IGF2BP2 and Promotes the Aerobic Glycolysis in Colorectal Cancer

The Role of lncRNA Crosstalk in Leading Cancer Metastasis of Head and Neck Squamous Cell Carcinoma

The Epitranscriptome beyond m6A

Potential Association between Asthma, Helicobacter pylori Infection, and Gastric Cancer

Over Expression of METRN Predicts Poor Clinical Prognosis in Colorectal Cancer

Prognostic and Predictive Value of m6A "Eraser

Regulation of N6-Methyladenosine in the Differentiation of Cancer Stem Cells and Their Fate

Systematic Screening of Protein-Coding Gene Expression Identified HMMR as a Potential Independent Indicator of Unfavorable Survival in Patients with Papillary Muscle-Invasive Bladder Cancer

Conversion of ATP to Adenosine by CD39 and CD73 in Multiple Myeloma Can Be Successfully Targeted Together with Adenosine Receptor A2A Blockade

Dynamic Transcriptomic m6A Decoration: Writers, Erasers, Readers and Functions in RNA Metabolism

Oncogenic Role of Long Noncoding RNAMALAT1 in Thyroid Cancer Progression through Regulation of the miR-204/IGF2BP2/m6A-MYC Signaling

Metformin Enhances Gefitinib Efficacy by Interfering with Interactions between Tumor-Frontiers in Cell and Developmental Biology | www

ILRC: a Hybrid Biomarker Discovery Algorithm Based on Improved L1 Regularization and Clustering in Microarray Data

The Fragile X Mental Retardation Protein Regulates Tumor Invasiveness-Related Pathways in Melanoma Cells

Reversible RNA Modification N 1 -methyladenosine (M 1 A) in mRNA and tRNA

Clinical Significance and Inflammatory Landscapes of a Novel Recurrence-Associated Immune Signature in Early-Stage Lung Adenocarcinoma

Prognostic Value of the Preoperative Fibrinogen-To-Albumin Ratio in Pancreatic Ductal Adenocarcinoma Patients Undergoing R0 Resection

Roles of N6-Methyladenosine (m6A) in Stem Cell Fate Decisions and Early Embryonic Development in Mammals

A Novel Prognostic mRNA/miRNA Signature for Esophageal Cancer and its Immune Landscape in Cancer Progression

Integrative Analysis of Multi-Omics Identified the Prognostic Biomarkers in

Analysis of Genetic Alteration Signatures and Prognostic Values of m6A Regulatory Genes in Head and Neck Squamous Cell Carcinoma

Mutant P53 in Cancer Progression and Targeted Therapies

A Novel Hypoxic Long Noncoding RNA KB-1980E6.3 Maintains Breast Cancer Stem Cell Stemness via Interacting with IGF2BP1 to Facilitate C-Myc mRNA Stability

Up-Regulation of SLC26A6 in Hepatocellular Carcinoma and its Diagnostic and Prognostic Significance

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.718974/ full#supplementary-material 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.