RNA transcriptome data of 404 tumor samples and 54 normal samples were obtained (https://portal.gdc.cancer.gov/), as well as clinical information on patients with OSC. The clinical information included age, sex, race, smoking, alcohol, pathologic stage, stage N, stage T, and tumor grade. To reduce statistical bias, we excluded patients with no survival information or a survival time of less than 3 months from further evaluation.

After reading and summarizing the published literature, 23 m7G regulatory factors were identified. The m7G regulatory factors include METTL1, WDR4, NSUN2, DCP2, DCPS, NUDT10, NUDT11, NUDT16, NUDT3, NUDT4, NUDT4B, AGO2, CYFIP1, EIF4E, EIF4E1B, EIF4E2, EIF4E3, GEMIN5, LARP1, NCBP1, NCBP2, NCBP3, EIF3D, EIF4A1, EIF4G3, IFIT5, LSM1, NCBP2 L, and SNUPN (Rong et al., 2022). A total of 16,773 lncRNAs were obtained in the OSCC data set of TCGA. We constructed a coexpression network of m7G methylation regulatory factors and lncRNAs by using the “igraph” package in the R program, and the conditional screening was (Pearson R)>0.3, p < 0.001. These lncRNAs are significantly related to the m7G methylation regulator (m7GR) and have been identified as m7G methylation-related lncRNAs.

Differentially expressed genes (DEGs) were detected between normal and tumor tissue samples by FDR 0.05 and |log2FC| > 2. According to the descending value of |log2FC|, the expression data of upregulated DEGs and downregulated DEGs were used, and R 4.1.2 was used to perform two-way hierarchical clustering analysis.

To identify lncRNAs with prognostic significance, univariate Cox regression analysis was performed on tumor samples and the control group. By using LASSO regression, a m7G-related lncRNA prognostic model (m7G-LncM) for OSCC patients was constructed using the R package glmnet. The steps for calculating the LASSO regression model’s risk score are as follows: Risksore = ∑ i = 1 n β i × Ε i

Ei is the expression value of the i gene in the model, and βi is the coefficient calculated by LASSO. According to this equation, we calculated the risk score of each OSCC patient. Patients were divided into high-risk and low-risk groups based on their median risk scores. Risk scores were compared to the overall survival rate (OS) at 1, 3 and 5 years using the receiver operating characteristic (ROC) curve, which was used to determine the area under the curve (AUC). To further demonstrate that the risk model of m7G-related lncRNAs has the ability to distinguish tumor patients, the “Limma” and “scatterplot3D” packages in R were used to conduct principal component analyses (PCA) on the genes associated with risk and the lncRNAs associated with m7G.

An independent prognostic factor included in the nomogram was screened using univariate and multivariate Cox regression (Wang et al., 2021). In accordance with the results of the univariate and multivariate Cox regressions, the probability of 1-, 3- and 5-year OS of OSC patients was evaluated with the “RMS” R package. The C-index and calibration curve were used to assess the nomogram’s capacity for prediction and discrimination. Decision curve analysis (DCA) was used to evaluate the net benefits of patients with OSCC under different threshold probabilities. The completion of DCA depends on two packages (stdca and dca).

For the score calculation of immune cell infiltration, ssGSEA was carried out with the gsva software package to calculate immune cell scores and assess immune-related pathways (p < 0.05) (Chen et al., 2022).

After obtaining the nucleotide variation data of OSCC, the genetic mutation burden of OSCC samples was analysed by Perl software. In OSCC patients with elevated and low risks for tumor mutation, the “limma” software package was used to analyze the difference in tumor mutation burden. Last but not least, all OSCC patients were grouped according to tumor mutation burdens and survival analysis was combined with the previously analysed high- and low-risk groups. With the survival and survminer R packages, we conducted survival analysis of mutation burden analysis.

To explain the differential signaling pathways between the high-risk group and the low-risk group of OSCC patients, we further understood the biological functions of these pathways. The method used in this study was Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

Based on the large-scale drug screening data in the database of cancer drug sensitivity genomics (GDSC) (https://cancerrxgene.org), this study combined with this genome analysis is beneficial to systematically identify the tumor response to chemotherapy drugs. In addition to calculating the maximum inhibitory concentration (IC50) of commonly used chemotherapy drugs for OSCC, an evaluation of its clinical utility was conducted. The R package used in this operation is “pRRophetic.” An analysis of Wilcoxon symbolic rank tests was carried out to compare the differences in IC50 between the high-risk and low-risk groups. For the data result, a block diagram is created by “pRRophetic” and “ggplot2” (Geeleher et al., 2014).

This study was conducted with tissues kindly provided by the Department of Stomatology of the Affiliated Hospital of Chengde Medical University. Medical ethics approval was obtained for sample acquisition. A total of three nontumor oral tissues and three OSCC tissues were obtained in October 2022 and stored at −80°C. Total RNA was isolated with TRIzol reagent, and quantitative real-time fluorescence quantitative PCR (qRT‒PCR) was performed with primers from Suzhou Genewiz Company and premix from SYBR Master Mix. The results were analysed using the 2−ΔΔCt method, and each sample test was repeated three times. Supplementary Table S1 contains details about the primer sequences.

First, according to Pearson correlation criteria, m7G-related lncRNAs were screened between 404 OSCC patients and 54 normal tissues. There were 1302 lncRNAs in total. In Figure 1A, the lncRNAs related to the m7G gene are shown as a network diagram. A total of 206 long noncoding RNAs differed in expression between the original cancerous tissue and normal tissue according to a differential gene expression analysis (Figure 1B), of which 36 genes were downregulated and 170 genes were upregulated (Figure 1C).

Univariate Cox regression was applied to 206 differentially expressed lncRNAs, and 11 m7G-related lncRNAs were identified (Figure 2A). To demonstrate this network of risk types related to the m7G regulator and lncRNAs, a Sankey diagram was constructed. (Figure 2B). Then, 11 m7G-related lncRNAs were subjected to Lasso regression to obtain five m7G-related lncRNAs. According to five m7G-related lncRNAs (Figures 2C,D), the prediction signal (figure) was constructed, and the risk coefficient of each lncRNA was obtained. The formula calculates the risk score as follows: AC108488.3*(−0.635)+AL133444.1*(0.683)+AC007128.1*(1.056)+AL359091.4* (0.382)+AL162413.1*(0.004). In addition, according to the model of all genes, m7G-related genes, m7G all-related lncRNAs and five m7G-related lncRNAs, PCA was conducted in this study. High-risk vs. low-risk differences among five m7G-related lncRNAs were the most obvious among the five (Figures 2E–H). As a result, the model we built is capable of distinguishing people at high and low risk.

According to the coefficient value of lncRNAs, it can be judged whether a lncRNA is a risk or a protection. Following the risk score, each OSCC sample from TCGA was divided into two low- and high-risk groups (Figure 3), using Kaplan-Meier survival curves to reflect the relationship between the risk score and survival rate. As shown in Figure 4A, there was a significant difference between the two groups (p < 0.001). The risk curve and graph indicate that the number of deaths increases obviously with increasing risk score (Figure 3A). This study showed a good correlation between OSCC survival and risk score. In addition, the multitime ROC curve demonstrated that risk signatures accurately predicted OSCC survival rates at 1, 3, and 5 years (Figures 4B–D). These results indicate that these five m7G-related lncRNAs play an important biological role in the development of OSCC. Additionally, this study compared the clinical characteristics of high- and low-risk patients (Figure 5). The results There was no difference in any kind of clinical information among the different groups, except the status condition.

As shown in Figure 6, considering the clinical information of OSCC patients (including age and sex) and m7G-related lncRNAs, univariate Cox regression (Figure 6A) and multivariate Cox regression (Figure 6B) were performed. Finally, three variables were included in this study as a nomogram to construct prognosis (Figure 6C). The calibration curve’s 1, 3, and 5-year predictions demonstrate that the development of the nomogram has some predictive value for the prognosis of OSCC (Figure 6E). In addition, OSCC prognosis is improved by the C index over other clinical information. In DCA, it can be seen that the nomogram has more clinical value for the prognosis of OSCC patients at 1, 3 and 5 years than other clinical stages (Figures 6F–H).

Regarding functional analysis, our study also examined 16 immune cells and 13 immune-related pathways for enrichment fraction and activities between low-risk and high-risk OSCC patients by single-sample genome enrichment analysis (ssGSEA). In patients with OSCC (Figure 7A), it was usually the high-risk subgroup that had high levels of immune cell infiltration, especially CD8⁺ T cells. Low-risk groups had downregulated APC coinhibition, cytolytic activity, and HLA response pathways. (Figure 7B).

According to the KEGG pathway analysis, the high-risk group did not have an enrichment in the signaling pathway, while chronic myeloid leukemia, dorsoventral axis formation, drug metabolism cytochrome p450 and the VEGF signaling pathway were enriched in the low-risk group (Figure 7C). This indicates that the signaling pathway and its related lncRNAs affect the occurrence, development and prognosis of OSCC, which is a mechanism worthy of further study.

Our first step in examining mutation differences between high-risk and low-risk groups was to examine mutation rates. The mutation frequency in the high-risk group (92.45%) was lower than that in the low-risk group (96.08%), according to the results (Figures 8A,B).In both groups, the most frequently mutated gene was TP53. In this study, we collected the somatic mutation data of OSCC and calculated the corresponding TMB score to study the potential role of OSCC tumor mutation load. It can be seen in Figure 8E shows that the mutation load in OSCC did not differ significantly between the high- and low-risk groups. With the median as the boundary point, patients were divided into “high TMB” and “low TMB” groups, and the results of survival analysis showed that there was no significant difference between the two groups (Figure 8C). It is possible, however, to determine the combined survival curve by combining tumor mutation load and risk score analysis. An interesting finding was that TMB and risk score were significantly correlated with OSCC survival (Figure 8D).

To better understand the individualized treatment of OSCC with drugs, this study used the drug analysis of the pRRophetic package to treat OSCC to obtain anticancer drugs sensitive to OSCC. A significant difference in doxorubicin sensitivity was seen between the high-risk and low-risk groups, as shown by Figure 9. This finding indicates that doxorubicin has a potential effect in the treatment of OSCC.

To further verify the expression levels of five m7G-related lncRNAs in nontumor oral tissues and OSCC tissues, the expression levels of the five lncRNAs in OSCC tissues all showed an upwards trend after RT‒qPCR (Figure 10). Further statistical analysis, however, indicated that only AL133444.1 and AL359091.4 showed significant differences between nontumor tissues and OSCC. The other three lncRNAs (AC108488.3, AC007128.1 and AL162413.1) did not show a statistically significant difference. This should be attributed to the small sample size in this study.