{'fraction_non_alphanumeric': 0.044925720422409166, 'fraction_numerical': 0.02447646321818507, 'mean_word_length': 5.1567493112947655, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of m6A methyltransferase-related genes predicts prognosis and immune infiltrates in head and neck squamous cell carcinoma

Background Head and neck squamous cell carcinoma (HNSCC) accounts for 90% of head and neck malignant tumors. As the early symptoms of HNSCC are not obvious, and it is prone to recurrence and metastasis, making the overall survival (OS) rate of patients very low. Existing studies have shown m6A methylation plays a crucial role in various cancers, but it is rarely studied in HNSCC. This study aimed to explore the expression of m6A methylation-related genes in HNSCC and its correlation with prognosis, and to explore its relationship with immune infiltration. Methods The gene expression data of HNSCC patient tumor samples (tumor =510) and adjacent normal tissue samples (normal =50) were extracted from The Cancer Genome Atlas (TCGA) database, and the expression characteristics of m6A regulatory factors were described. Kaplan-Meier survival analysis was used to analyze the relationship between m6A regulatory factors and OS and disease-specific survival (DSS). Least absolute shrinkage and selection operator (LASSO) regression was used to construct the m6A regulatory factor-HNSCC risk prediction model. In addition, the relationship between m6A methylation-related genes and tumor immune infiltration were discussed. Results The differential expression of 20 genes were identified by TCGA, and 18 genes (IGF2BP2, IGF2BP1, IGF2BP3, VIRMA, YTHDF1, YTHDF2, YTHDF3, ZC3H13, METTL14, ALKBH5, METTL3, RBMX, WTAP, YTHDC1, FTO, HNRNPC, HNRNPA2B1, and RBM15) were overexpressed in HNSCC. The survival rate of different gene expression levels was different. The high expression of YTHDC1 and YTHDC2 indicated better OS. Furthermore, for DSS, increased expression of YTHDC2 was also correlated with better clinical outcomes (P<0.05). At the same time, we drew a 3-gene risk score model in the TCGA-HNSCC cohort, and the survival curve showed compared with low-risk patients, high-risk patients had significantly worse OS (P<0.05). Gene enrichment analysis showed EPITHELIAL_MESENCHYMAL_TRANSITIO, MTORC1_SIGNALING, MYC_TARGETS_V1, MYC_TARGETS_V2, MYOGENESIS pathways, high TP53 mutations, and suppressive immunity were related to the high-risk group. The low-risk group was related to ALLOGRAFT_REJECTION, COMPLEMENT, IL6_JAK_STAT3_SIGNALING, INTERFERON_ALPHA_RESPONSE, INTERFERON_GAMMA_RESPONSE pathways, low TP53 mutations, and active immunity. Conclusions The m6A methyltransferase-related genes can predict the prognosis of HNSCC and are related to immune infiltration.


Introduction
Head and neck squamous cell carcinoma (HNSCC) is a major global public health problem, with a high rate of recurrence and distant metastasis (1). It is a heterogeneous disease even though the tumors all originate in the epithelial cells of the mucosal linings of the upper airway and food passages (2). The onset of HNSC is hidden and the early clinical symptoms are not obvious, so patients are mostly in the middle and late stages of disease at the time of diagnosis. Studies have shown that the occurrence and metastasis of HNSCC is complex, involving multiple genes. Therefore, searching for genetic and molecular markers related to the occurrence and development of HNSCC will help to better evaluate its prognosis and find new treatment methods, thereby improving the survival rate of HNSCC patients. Meanwhile, it is very important to understand the tumor immune microenvironment of HNSCC and identify valuable biomarkers. One research has developed an immune-related genetic prognostic index (IRGPI), and after analyzing different molecular and immune characteristics. It is found that in the classification of HNSCC immune subtypes, IRGPI-high subtypes are higher than IRGPI-low subtypes. In the high subgroup of IRGPI, cytotoxic CD8 T cells, CD4 T cells and M1 macrophages are more enriched. B cells and M0 and M2 macrophages are more enriched in the IRGPI low subgroup.
The most common post-transcriptional modification in eukaryotic RNA is N6-methyladenosine (m6A) (3). This modification mainly occurs in the DRACH (D = A, G, U; R = A, G; H = A, C, U) in various RNAs with adenine in the sequence (4). As early as the 1970s, scientists discovered the m6A modification, which is significantly concentrated in the messenger RNA (mRNA) stop codon and 3' noncoding region (5), which can dynamically regulate biological processes such as cell differentiation. Although m6A modification has been discovered for more than 40 years, it is difficult to detect due to its poor stability and short life span (6). The current research found that m6A is a reversible dynamic process, and its function is jointly determined by the "Writer", "Eraser", and "Reader". The "Writer" is also called methyltransferase complex, the known components of which include METL3 (7), METL14, and WTAP (8,9). The "Erase", also called demethylase, including ALKBH5 (10) and FTO (11), can reverse methylation n. The functions of m6A are recognized by the m6A binding protein, which is the "reader". The currently discovered m6A binding proteins include the YTH (YT521-B homology) domain protein family, which is YTHDF1, YTHDF2, YTHDF3 (12), YTHDC1 (13), and YTHDC2 (5) and the HNRNP family of heterogeneous nuclear proteins, namely HNRNPA2B1 and HNRNPC (14).
The evolutionarily conserved YTH domain contained in the YTH protein family can selectively recognize m6A (13). The family member YTHDF1 promotes the binding of target RNA to ribosomes by interacting with initiation factors, and enhancing mRNA translation and protein synthesis (15). Additionally, YTHDF2 reduces mRNA stability and accelerates mRNA degradation by recruiting CCR4-NOT complexes (12,16). The effect of YTHDF3 is bidirectional. It can not only interact with YTHDF1 to enhance RNA translation but also combine with YTHDF2 to promote RNA degradation (17). A function of YTHDC1 is to affect the splicing and output of RNA, and YTHDC2 can promote the translation efficiency of target RNA (18).
More and more evidences show that m6A mRNA methylation is closely related to the occurrence and development of tumors, and m6A modification has different biological functions in different types of tumors. Glioblastoma stem cells (GSCs) are a kind of stem cells with strong invasive ability. In recent years, studies have found that GBM is closely related to the level of m6A increases and the induced differentiated of m6A moCs. Inhibiting the expression of METTL3 or METTL14 can reduce m6A levels, enhance the growth and self-renewal of GSCs in vitro, and promote the ability of GSCs to form tumors in vivo. In addition, METL3 is significantly up-regulated in human hepatocellular carcinoma cells. Overexpression of METL3 can promote HCC proliferation and migration, and significantly enhance HCC tumor formation. Suppressor of cytokine signaling 2 (SOCS2) is a downstream target of METTL3 and also a tumor suppressor. METTL3 reduces the stability of SOCS2 mRNA through the m6A-YTHDF2-dependent pathway to promote the occurrence and development of HCC. Studies have shown that m6A methylation-related genes are involved in the occurrence and development of many diseases, but there is very little research regarding their role in in HNSCC. The purpose of this study was to explore the effect of m6A methylation-related genes on the prognosis of HNSCC and their relationships with tumor immune infiltration. We present the following article in accordance with the REMARK reporting checklist (available at https:// dx.doi.org/10.21037/atm-21-4712).

The TCGA-HNSCC cohort data collection
The RNA-seq data (Fragments Per Kilobase per Million (FPKM)) of the HNSCC dataset, corresponding clinical information, and somatic mutation data were obtained from TCGA GDC (https://portal.gdc.cancer.gov/). The RNA sequencing data (FPKM values) were then transformed to transcripts per million reads (TPM) and normalized into log2 (TPM +1). Subsequently, data analysis was performed in R language (R version 3.6.3; https://cran.r-project.org/bin/windows/base/old/3.6.3/). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Tumor-infiltrating immune cells corresponding to TCGA-HNSCC cohort
The Immune Cell Abundance Identifier (ImmuCellAI) is an online tool to estimate the abundance of 24 immune cells from gene expression dataset, in which the 24 immune cells are comprised of 18 T-cell subtypes and 6 other immune cells: B cell, natural killer (NK) cell, monocyte cell, macrophage cell, neutrophil cell, and dendritic cell (DC) (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) (19). The infiltrating data of tumor-infiltrating immune cells (TIICs) corresponding to TCGA-HNSCC samples were downloaded from the ImmuCellAI website.

Construction of the m6A methylation-related gene signature
The least absolute shrinkage and selection operator (LASSO)-Cox regression model was performed for variable screening and complexity adjustment to screen potential genes to establish the m6A methylation-related gene signature as previously described. The penalty value parameter was determined after 10-fold cross-validation using the "glmnet" package (20). Finally, a formula was developed using the gene signature constructed above. Risk score = expression of a gene (YTHDC2) × corresponding coefficient (−0.297) + expression of a gene (HNRNPA2B1) × corresponding coefficient (−0.752) + expression of gene (RBM15B) × corresponding coefficient (−0.332) (21,22). Patients of the TCGA-HNSCC cohort were divided into high and low-risk score groups.

Gene set enrichment analysis
To assess the potential differences in biological functions between the high-and low-risk score subgroups, gene set enrichment analysis (GSEA) software (https://www. gsea-msigdb.org/gsea/login.jsp) was used based on the hallmarks gene set ("h.all.v7.0.symbols.gmt") as previously described (22).

Additional bioinformatic and statistical analysis
We used R software (version 3.6.3, https://www.r-project. org/) to perform analyses and plot graphs. Wilcoxon's rank-sum test was used to analyze the parameters of different groups. Correlations between the levels of m6A methylation-related genes in HNSCC were also analyzed using Spearman correlation. A consensus cluster consisting of 20 m6A methylation-related genes was constructed by using the "ConsensusClusterPlus" package. The ImmuneScore was developed by using the "ESTIMATE" packages. The "survival" package was used to perform survival analysis (https://CRAN.R-project.org/ package=survival). The best critical value was determined using the R packages ("survival" and "survminer") and a double-sided log-rank test. Kaplan-Meier survival curves and log-rank tests were used to compare the differences in survival between high-and low-risk score groups. The R package "survivalROC" was used to plot the timedependent receiver operating characteristic (ROC) curve to determine the prognostic value of the signature with OS in HNSCC. Visualization of analysis results was performed using the R packages ggplot2 (http://ggplot2.org). A P value of <0.05 was considered statistically significant.

Patient characteristics
We downloaded RNAseq and clinical data in HNSCC tumor and normal tissues from the TCGA dataset. After excluding cases with insufficient OS time, or lack of data, our study included a total of 502 patients (368 males and 134 females). The median age of the participants was 61 years. Clinical data were collected including age, gender, grade, race, clinical stage, history of smoking and alcohol, tumor, node, metastasis (TNM) stage, and pathological grade (Table S1).

The ImmuneScore in differentially expressed m6A methylation-related genes level in HNSCC
The ImmuneScores were calculated through the ESTIMATE evaluation method, and the immune score in differentially expressed m6A methylation-related gene groups were analyzed (Figure 2, Figure S1). As shown in Figure 2A, patients with high expression of YTHDC1 had lower immune scores than those with low YTHDC1. In addition, 9 m6A methylation-related genes were observed to be negatively correlated with immune score, including IGF2BP1, IGF2BP3, YTHDF1, VIRMA, METTL3, HNRNPA2B1, ALKBH5, IGF2BP2, and HNRNPC (P<0.05). Immune score showed no differences in 8 differentially expressed m6A methylationrelated genes groups ( Figure S1) (P>0.05).

The prognostic values of m6A methylation-related genes in HNSCC
We next performed overall survival (OS) analysis for m6A methylation-related genes in HNSCC. The HNSCC patients with higher expression of HNRNPA2B1, HNRNPC, METTL3, IGFBP2, IGFBP1, and IGFBP3 had a poorer prognosis; however, those with a higher expression of YTHDC1 and YTHDC2 had a better prognosis ( Figure 3). Considering the possibility of non-tumor death during follow-up, we analyzed the relationship between m6A methylation-related genes expression and disease-specific survival (DSS) in HNSCC. For DSS, increased expression of YTHDC2 was correlated with better clinical outcomes in HNSCC patients ( Figure 4). Increased expression of IGFBP2, METTL3, HNRNPA2B1, and HNRNPC genes were correlated with poor clinical outcomes in HNSCC patients. Taken together, these data suggested that partially m6A methylation-related genes could be used as a biomarker predicting unfavorable prognosis in HNSCC patients.

Construction of the 3-m6A methylation-related genes-based classifier
We used a LASSO Cox regression model to build a superior prognostic classifier, which selected 3 m6A methylationrelated genes from 4 m6A methylation-related genes (IGF2BP2, YTHDC2, HNRNPA2B1, and RBM15B) identified in the TCGA-HNSCC cohort: YTHDC2, HNRNPA2B1, and RBM15B ( Figure 6A,6B). Using the LASSO Cox regression models, we calculated a risk score for each patient based on their individual expression levels of the 3 m6A methylation-related genes, where risk score = (−0.297 × status of YTHDC2) + (−0.752 × status of HNRNPA2B1) + (−0.332 × status of RBM15B). As shown in Figure 6C, compared with low-risk patients, high-risk patients had significantly worse OS (P<0.05). The timedependent ROC curve showed that the area under the curve (AUC) of the 3-, 5-, and 8-year OS in the cohort were 0.62, 0.603, and 0.643, respectively ( Figure 6D). Thus, the 3-m6A methylation-related genes-based classifier performed well in predicting OS in HNSCC patients.

Immune characteristics of 3-m6A methylation-related genes-based classifier subgroups
Wilcoxon test was used to analyze the composition of immune cells in 3-m6A methylation-related genes-based classifier subgroups. We found that B cells, CD4 T cells, CD8 T cells, iTreg, NK, Tc, Tcm, Tex, Tfh, Tgd, Th1 cells, Th2 cells, and Tr1 cells were more abundant in the lowrisk score subgroup, while CD8 naive cells, CD4 naive cells, DC, monocytes, neutrophils, natural killer T cells (NKT), and Th17 cells were more abundant in the high risk-score subgroup (Figure 7).

Molecular characteristics of patients in the high-and low-risk score subgroup
We used GSEA software to find the difference between the high-risk score subgroup and the low-risk score subgroup in the Hallmark pathway. As shown in Figure 8A Figure 8C,8D). Among them, the mutation rates of TP53, TIN, and FAT were greater than 20% in both the high-risk group and the low-risk group of HNSCC patients. In addition, TP53 mutations were the

Discussion
The occurrence and metastasis of HNSCC is a complex process involving multiple genes. Therefore, searching for genes and molecular markers related to HNSCC will help us better evaluate the prognosis of HNSCC patients. There is only one study that has reported a correlation between m6A and HNSCC, in which METTL3-mediated inhibition of ZNF750 promoted the progression of NPC. Previous studies have shown that METTL3 can promote apoptosis of a variety of tumor cells, such as breast cancer and gastric cancer. Therefore, in order to predict the prognosis of HNSCC patients, we used TCGA to establish 20 gene characteristics. The expression of all genes was positively correlated, and 18 m6A methylation-related genes were significantly up-regulated in HNSCC (including IGF2BP2,  IGF2BP1, IGF2BP3, VIRMA, YTHDF1, YTHDF2,  YTHDF3, ZC3H13, METTL14, ALKBH5, METTL3, RBMX, WTAP, YTHDC1, FTO, HNRNPC, HNRNPA2B1, and RBM15) (P<0.05). In addition, we also analyzed the relationship between these 18 m6A methylation-related genes and immune scores. The study showed that 10 m6A methylation-related genes were negatively correlated with immune scores (including YTHDC1, IGF2BP1, IGF2BP3, YTHDF1, VIRMA, METTL3, HNRNPA2B1, ALKBH5, IGF2BP2, and HNRNPC) (P<0.05). Next, the Gene Expression Profiling Interactive Analysis (GEPIA) database was used to analyze the most differential survival genes of HNSCC. The differential expression of the remaining eight m6A methylation-related genes did not differ from the immune score (P>0.05). Our research also showed that YTHDC2 and YTHDC1 can be used as protective genes because their high expression suggests a better prognosis, while the IGFBP2, METTL3, HNRNPA2B1, and HNRNPC genes are associated with and indicative of a poor prognosis.
Our results show that the expression levels of YTHDC1 and YTHDC2 are positively correlated with the prognosis of HNSCC and are expected to be prognostic biomarkers of HNSCC. These findings indicate m6A modulators can According to the expression levels of 20 m6A methylationrelated genes, we used R's ConsensusClusterPlus package to group the data into 3 subgroups. The overall analysis showed the highest and lowest risk ratios of subgroup 1 and subgroup 2 stages, indicating that the dynamic time is related to the comprehensive expression level of m6A methylation-related genes.
The LASSO algorithm simultaneously analyzes all independent variables and selects the most influential variable (21). According to LASSO Cox analysis, 3 m6A methylation-related genes (YTHDC2, HNRNPA2B1, and RBM15B) were selected from the 4 m6A methylationrelated genes (IGF2BP2, YTHDC2, HNRNPA2B1, and The results indicate that m6A methylation-related genes can independently predict the prognosis of HNSCC patients. Studies have shown that YTHDC2 can specifically recognize and bind m6A through the YTH sequence, promote mRNA translation and affect mRNA stability (4). The YTHDC2 gene is essential for the maintenance of meiosis and the production of germ cells (23)(24)(25), and is involved in the occurrence and development of a variety of diseases, including autism, pancreatic cancer, colon adenocarcinoma, liver cancer, and lung adenocarcinoma (26)(27)(28)(29)(30). The HNRNPA2B1 gene was highly expressed in a variety of human cancers, such as prostate cancer (31), pancreatic cancer (32), and hepatocellular carcinoma (33). It has been reported that RBM15B participates in the evolution of a variety of solid tumors, such as when Zhang et al. (34) found that RBM15B has at least 2 mutated possible carcinogenic sites in primary epithelial ovarian cancer through next-generation highthroughput sequencing. Shahriyari et al. (35) found that it was highly positively correlated with the expression of RBM15B. The BRCA1 associated protein-1 (BPA1) is a prognostic gene for uveal melanoma and breast cancer. In addition, RBM15B is also associated with an increased risk of malignant mesothelioma and renal clear cell carcinoma (26,36,37). This is consistent with our research results.
Understanding the status of the tumor microenvironment (TME) can help find new ways to treat HNSCC or change TME to improve the effectiveness of immunotherapy. We divided the 3 m6A methylation-related genes into high-and low-risk groups and analyzed the composition of immune cells between the 2 subgroups. The results showed that B cells, cytotoxic CD8 T cells, and CD4 T cells were more abundant in the low-risk subgroup, while inflammatory cells such as neutrophils, monocytes, and DC cells were more common in the high-risk subgroup. A large number of studies have shown that dense infiltration of T cells, especially cytotoxic CD8 T cells, indicates a good prognosis (38)(39)(40). While tumors stimulate neutrophils to promote tumor migration, invasion, and metastasis (41), acute and chronic inflammation are conducive to tumor growth and the development of aggressive phenotypes. Our research results support these conclusions. In addition, GSEA showed that the high-risk subgroup of m6A methylation-related genes has immunosuppressive effects and is related to high TP53 mutations, while the lowrisk subgroup is related to active immunity and low TP53 mutations.
Our research has several limitations. First of all, the sample size of this study is small, and it only focuses on the mRNA sequencing data from TCGA. Future work is necessary to further improve the sample size of HNSCC patients and further explore other public databases to obtain more verification from independent cohorts to confirm the results. Secondly, our conclusion lacks certain clinical practice, so more clinical studies are needed to confirm whether m6A-related genes can be used as new targets for HNSCC treatment. In conclusion, our study has determined the risk characteristics of m6A methylation-related genes, which can be independently used to predict the prognosis of HNSCC patients and is related to immune infiltration, providing a new treatment strategy for HNSCC patients.  **** **** *** **** **** **** **** **** **** **** **** **** **** **** ** ** **   Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the noncommercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.