m6A-Related lncRNAs Predict Overall Survival of Patients and Regulate the Tumor Immune Microenvironment in Osteosarcoma

Background m6A-related lncRNAs have demonstrated great potential tumor diagnostic and therapeutic targets. The goal of this work was to find m6A-regulated lncRNAs in osteosarcoma patients. Method The Cancer Genome Atlas (TCGA) database was used to retrieve RNA sequencing and medical information from osteosarcoma sufferers. The Pearson's correlation test was used to identify the m6A-related lncRNAs. A risk model was built using univariate and multivariable Cox regression analysis. Kaplan–Meier survival analysis and receiver functional requirements were used to assess the risk model's performance (ROC). By using the CIBERSORT method, the associations between the relative risks and different immune cell infiltration were investigated. Lastly, the bioactivities of high-risk and low-risk subgroups were investigated using Gene Set Enrichment Analysis (GSEA). Result A total of 531 m6A-related lncRNAs were obtained from TCGA. Seven lncRNAs have demonstrated prognostic values. A total of 88 OS patients were separated into cluster 1, cluster 2, and cluster 3. The overall survival rate of OS patients in cluster 3 was more favorable than that of those in cluster 1 and cluster 2. The average Stromal score was much higher in cluster 1 than in cluster 2 and cluster 3 (P < 0.05). The expression levels of lncRNAs used in the construction of the risk prediction model in the high-risk group were generally lower than those in the low-risk group. Analysis of patient survival indicated that the survival of the low-risk group was higher than that of the high-risk group (P < 0.0001) and the area under the curve (AUC) of the ROC curve was 0.719. Using the CIBERSORT algorithm, the results revealed that Macrophages M0, Macrophages M2, and T cells CD4 memory resting accounted for a large proportion of immune cell infiltration. By GSEA analysis, our results implied that the high-risk group was mainly involved in unfolded protein response, DNA repair signaling, and epithelial-mesenchymal transition signaling pathway and glycolysis pathway; meanwhile, the low-risk group was mainly involved in estrogen response early and KRAS signaling pathway. Conclusion Our investigation showed that m6A-related lncRNAs remained tightly connected to the immunological microenvironment of osteosarcoma tumors, potentially influencing carcinogenesis and development. The immune microenvironment and immune-related biochemical pathways can be changed by regulating the transcription of M6A modulators or lncRNAs. In addition, we looked for risk-related signaling of m6A-related lncRNAs in osteosarcomas and built and validated the risk prediction system. The findings of our current analysis will facilitate the assessment of outcomes and the development of immunotherapies for sufferers of osteosarcomas.


Introduction
Osteosarcoma (also known as osteogenic sarcoma or simply bone cancer) is the most frequent kind of bone cancer [1]. It is a type of invasive malignant neoplasm that develops from primitive altered cells of mesenchymal origin (and hence is a sarcoma) and demonstrates osteoblastic differentiating and creates cancerous osteoid [1]. e most frequent histological kind of primary bone sarcoma is osteosarcoma [2]. It is particularly common among adolescents and young adults [3].
Despite the positive results of comprehensive surgical resection along with chemotherapy and radiotherapy, roughly 40-50 percent of patients develop lung dissemination [4]. Patients with lung metastatic tumors had a five-year life expectancy of only 28% [5]. erefore, searching for new treatment targets and prognostic biomarkers is critical.
N6-methyladenosine (m6A) is the methylated adenosine that occurs at the N6 site and it is prevalently an epigenetic alteration in mRNA and noncoding RNAs (ncRNAs) [6]. e process of m6A is considered to be dynamic and reversible. Similar to other epigenetic regulatory mechanisms, the bioactivities of m6A are coordinated by "writer," "reader," and "eraser" [7]. e compound of m6A-writer acts as a methylase, allowing m6A to be installed. e m6A alteration is then identified by m6A-associated proteins, which are also called as readers. e erasers are demethylases that are in charge of eliminating the m6A alterations. Several investigations have shown that m6A alterations control carcinogenesis and progression [8]. e relevance of m6A variations in malignancy prediction is becoming clear.
Long noncoding RNA (lncRNA) is RNA with a transcriptome sequence of more than 200 nucleotides that are not transcribed into proteins [9]. lncRNA is involved in a number of physiological systems in cells, including cell growth and differentiation [10]. lncRNA influences gene transcriptions at three different sides: epigenetic regulation, transcription, and posttranscription. e abnormal pattern of lncRNA is linked to cancerous aggressiveness as well. Recent evidence has proved that lncRNA acts as tumorigenesis or repressors by participating in numerous signaling pathways [11]. For instance, SRY-box transcription factor 2 (SOX2) promotes development in colorectal cancer (CRC) by catalyzing methyl with methyltransferase-like protein 3 (METTL3) [12], whereas BCL2 interacting protein 3 (BNIP3) promotes cancer progression in breast cancer by catalyzing demethylation with FTO (fat mass and obesityassociated protein) [13]. Furthermore, lncRNAs have the potential to be used as prognosis markers in various types of tumors, such as prostate, breast, and liver cancers [14]. Recently, there has been a growing interest in the research framework between m6A and lncRNAs. eir regulating complex is implicated in tumor growth, migration, and metastasis in a variety of malignancies, giving new targets for cancer diagnostics and treatment. lncRNAs and m6A have been shown to have critical roles in controlling cancer bioactivity, but the mutual regulation mechanism between them is unknown. As a result, it is critical to investigate the possible biochemical reaction of these unregulated m6Arelated lncRNAs in osteosarcoma. erefore, determining the associations between the alterations of m6A-related lncRNAs and the courses of osteosarcoma may help to find biomarkers that can be considered predictive and prognostic markers.

Screening of Prognostic Genes.
In accordance with the TCGA naming convention, all samples ending in 01A were OS samples, including 88 cases in total in our study. Gene IDs for TCGA-OS expression profiles were converted using the package (org.Hs.eg.db), obtaining a total of 34446 gene symbols. When there was a single symbol corresponding to multiple IDs, the expression profiles of identical gene symbols were merged by the maximum value. e GRCh38 annotation file was downloaded from the GENCODE website, which was used to differentiate the types of genes for TCGA-STAD expression profiles. e expression profiles of 20 m6A-related genes were extracted and box-line plots were plotted using the ggplot2 package. Correlation analysis was conducted on the expression profiles of all samples using the rcorr function of the Hmisc package. e correlation cor values and p values of the 10 m6A-related genes with lncRNA genes were obtained using Pearson analysis (p value < 0.05, cor Filter>0.2) [16]. A total of 531 m6Arelated lncRNAs were obtained.

Cox Regression Analysis.
Clinical data packages for survival time and survival status were combined with 531 lncRNA gene expression profiles associated with m6A. Oneway Cox regression analysis was performed using the R package survivor (p value < 0.05). Correlation analysis was carried out using the cor.test of the stats package (p value < 0.05).

Consistent Clustering Analysis.
e expression profiles of seven lncRNA prognostic genes associated with disease risk were analyzed by consistent clustering using the R package consensus cluster plus, setting the number of clusters to 2.

Prognosis and the Immune Microenvironment in Coherent
Clustering.
e expression profiles of 9 prognostic genes were mapped using the pheatmap package, plus additional grouping tags (such as age and gender), and divided into Cluster 1, Cluster 2, and Cluster 3. Survival curve analysis was performed with the survfit of the surv. package and plotted with the ggsurvplot function of the survminer package, grouped as Cluster1, Cluster2, and Cluster3. Using the CIBERSORT algorithm, immune infiltration of patient tissue is identified by 22 different types of immune cells. LM22 feature matrix file with CIBERSORT algorithm (1000 permutations) was applied to compare the immunocyte infiltration scores of Cluster1, Cluster2, and Cluster3.
Using the ESTIMATE package, immune scoring was calculated for TCGA-OS expression profiles, including ESTI-MATE score, Immune score, Stromal score, and Tumor Purity. e box plots were drawn using the ggplot2 package, grouped as Cluster1, Cluster2, and Cluster3. e significance of differences between groups was calculated by T test.

Constructing the Prognosis Diagnosis Model of TCGA-OS.
Risk scores were calculated for the seven lncRNAs screened that were associated with OS risks. e calculation formula is as follows: where coef(lncRNAn) represents the regression coefficient of lncRNA. Based on this formula, the risk score of each group could be derived [17]. e median risk score was used to classify the OS samples containing survival information into a high-risk group and a low-risk group. e risk score distribution of the sample was plotted using the ggplot2 package. e samples were divided into high-risk and low-risk groups. e survival curves were analyzed using the survfit function of the survivor package and plotted using the ggsurvplot function of the survminer package.

Immune Infiltration Analysis.
e LM22 feature matrix file with the CIBERSORT algorithm (1000 permutations) was used to compare the immunocyte infiltration scores of the high-risk group and low-risk group. e risk core risk scores of Cluster1, Cluster2, and Cluster3 were compared, and differences between groups were calculated using a ttest.
e box plot was obtained by using the ggplot2 package.

Correlation of Immune Infiltration with lncRNA.
e correlation between hub gene expression and immune infiltrating cells was calculated by the cor.test function of the psych package (spearman algorithm). Heatmaps were drawn with the pheatmap package to show the correlation heatmaps ( * * * P < 0.001, * * P < 0.01, and * P < 0.05). e colors from blue to red indicate a correlation cor-value from small to large.

e m6A-Related Differential Genes and Identification of m6A-Related lncRNAs in OS.
88 samples of OS were derived from the TCGA database. e samples were divided into two subgroups: metastases (22 specimens) and nonmetastases (65 specimens). 3039 lncRNAs and 19135 mRNAs were obtained in total. 17,249 differential genes were found using the Limma package differential analysis (p value < 0.05).
28 lncRNA prognostic genes associated with the risk of OS were found in total. Multifactor Cox regression analysis was continued for 28 lncRNAs obtained by single-factor regression using the R package SURVIVAL (p value < 0.05) [17]. A total of 7 lncRNA prognostic genes were screened for association with prognostic significance (Table 1). e expression profiles of 10 differential genes regulated by m6A and 7 lncRNA prognostic genes associated with risk factors were extracted separately. e correlation between m6Arelated genes and lncRNAs is shown in Figure 1.

Consensus Clustering Categorized OS Patients according to m6A-Related Prognostic lncRNAs.
Following the transcription of m6A-related predictive lncRNAs, consensus clustering was also used to divide OS patients into different clusters ( Figure 2). Due to the similarities demonstrated by the expression profiles of m6A-related predictive lncRNAs, k � 3 was discovered to be the best clustering consistency from k � 2 to 9 (Figure 2(a)). A total of 88 OS patients were separated into cluster 1, cluster 2, and cluster 3. e overall survival rate of OS patients in cluster 3 was more favorable than that of those in cluster 1 and cluster 2 ( Figure 2

Consensus Clustering Linked to Immune Infiltration.
To explore the role of m6A-related prognostic lncRNAs in the osteosarcoma immune microenvironment, we then analyzed the difference in the immune score and immune cell infiltration level among cluster 1, cluster 2, and cluster 3. e average Stromal score was much higher in cluster 1 than those in cluster 2 and cluster 3 ( Figure 3) (P < 0.05).

Risk Signaling of m6A-Related lncRNAs and Risk Prediction Model Showed Prognostic Value in Osteosarcoma.
After calculating the risk scores of individual patients, the cases were divided into a high-risk group and a low-risk group based on the median risk score (Figure 4(a)). e heatmap showed changes in the risk score and expression levels of lncRNAs. e expression levels of lncRNAs used in the construction of the risk prediction model in the high-risk group were generally lower than those in the low-risk group (Figure 4(e)). Analysis of patient survival indicated that the survival of the low-risk group was higher than that of the Computational Intelligence and Neuroscience high-risk group (P < 0.0001, Figure 4(d)). e area under the curve (AUC) of the ROC curve was 0.719, indicating that the model showed medium accuracy in predicting the prognosis of osteosarcoma patients (Figure 4(e)). Our results showed that m6A-related lncRNAs had prognostic value and that the risk prediction model showed satisfactory performance in predicting the prognosis of patients with osteosarcoma ( Figure 4).

Correlations between Immune Infiltration and m6A-Related lncRNAs in Osteosarcoma.
e infiltration of 22 immune cell types in osteosarcoma samples was further analyzed using the CIBERSORT algorithm. e results revealed that Macrophages M0, Macrophages M2, and Tcells CD4 memory resting accounted for a large proportion of immune cell infiltration ( Figure 5). However, there were no significant differences in the proportion of most of the infiltrated immune cell types between the high-and low-risk groups. e correlation coefficients between 22 types of immune cells are depicted in Table 2 and Figure 6. 3.6. e m6A-Related lncRNAs Subgroups Were Associated with Biological Characteristics of Osteosarcoma. Next, we analyzed differences in the biological responses of the subgroups generated using consensus clustering in order to further explore the relationships between m6A-related lncRNAs clustering subgroups and osteosarcoma biological characteristics. GSEA was used to explore the main KEGG signaling pathways in the two subgroups. Our results showed that the high-risk group was mainly involved in unfolded protein response, DNA repair signaling, and epithelial-mesenchymal transition signaling pathway and glycolysis pathway (Figures 7(a)-7(d)), which are related to the initiation of metastasis in cancer progression [18][19][20]. e low-risk group was mainly involved in estrogen response early and KRAS signaling pathway (Figures 7(e) and   Computational Intelligence and Neuroscience 7(f )), which are important for oncogene of osteosarcoma [21,22]. ese results indicated that there were differences in biological characteristics between high-risk and low-risk groups, which affected tumorigenesis and progression of osteosarcoma and survival of osteosarcoma patients.

Discussion
With the advancement of high-throughput sequencing methods and their widespread use in oncology, disclosing epigenomic irregularities during tumor onset and progression opens the door to the proof of identity of targeted therapies and predictive biomarkers in a range of tumors [23]. Current attention has been focused on the comprehensive detection of high sequencing results using publicly available databases, which has resulted in several major findings in the diagnosis and therapy of malignancies [17,24]. Conversely, the involvement of m6A regulators in lncRNA deregulation in human cancer is unknown. Few research studies investigated the associations between m6A alterations and lncRNA-dependent osteosarcoma.
To elucidate the underlying biological mechanisms of the m6A-related lncRNA profile, we identified genes that were expressed differently between the metastases and nonmetastases groups. A total of 7 lncRNA prognostic genes were screened for association with prognostic significance in osteosarcoma patients (including TNS1-AS1, WWC2-AS1, TFPI2-DT, LINC01474, LINC00910, LINC01982, and LINC00538). We did a clustering analysis based on the transcriptome of m6A-related predictive lncRNAs to further investigate the association between m6A-related lncRNAs and the clinical and pathological and biological aspects of AML. is analysis produced three subgroups, which were labeled as cluster 1, cluster 2, and cluster 3. e overall survival rate of OS patients in cluster 3 was more favorable than that of those in cluster 1 and cluster 2.
In addition, the samples were separated into two groups according to the average risk score: high-risk and low-risk groups after computing the risk scores of patient characteristics. GSEA suggested that the high-risk group was mostly engaged in unfolded protein response, DNA repair signaling, epithelial-mesenchymal transition signaling, and glycolysis signaling pathways. e low-risk group was      mainly involved in estrogen response early and KRAS signaling pathway. e low-risk group had a better overall survival rate than the high-risk group did. e ROC curve's area under the curve (AUC) was 0.719, meaning that the model was somewhat capable of predicting the prognosis of osteosarcoma patients. According to a prior study, the active  unfolded protein response was a key biochemical characteristic of osteosarcoma [25]. A large amount of evidence about the involvement of DNA damage response in osteosarcoma formation, survival, and treatments was assembled in a review paper [26]. Diverse immune cells and released substances make up the immunological microenvironment [27]. Tumor cell invasion, immunological capabilities, and the activation of targeted therapies could all determine cancer patients' outcomes and forecast their responsiveness to immunotherapies [28,29]. e immunological microenvironment and immune-related biochemical pathways could be altered by modulating the transcription of m6A regulatory or lncRNAs [30]. Our results revealed that Macrophages M0, Macrophages M2, and T cells CD4 memory resting accounted for a large proportion of immune cell infiltration. It is reported that M2 macrophages can encourage tumor progression [31,32]. Especially, elevated infiltration of M2 macrophages has been linked to osteosarcoma metastases and poor patient outcomes, notwithstanding the adoption of intensive therapy regimens [33]. Moreover, exosomes from metastasis osteosarcoma cells have been shown to alter tumor-associated macrophage cellular signaling, boost the M2 phenotype, and establish an immunosuppressive, tumor-promoting milieu via producing TGFB2 [34]. However, there were no significant differences in the proportion of most of the infiltrated immune cell types between the highand low-risk groups.
We investigated the association between m6A-related lncRNAs expression and the levels of immune infiltration in OS using CIBERSORT. TNS1-AS1 and TFPI2-DT were found positively correlated with B cell memory and B cells naive in our data, respectively. ese results about immune infiltrate in osteosarcomas were consistent with previous studies [35,36]. Additionally, the results uncovered that LINC01474 had a positive correlation with T cells CD8; however, LINC00910 was negatively related to T cells CD8. Other scholars believed that the significant infiltration of T cells CD8 activation in tumors may aid our signature's capacity to attain consistent predictive value [37]. ere has been proof that programmed death-1 receptor (PD-1) plays a role in the evolution of osteosarcomas, and the proportion of PD-1 in blood CD8+ T cells is much higher in osteosarcomas patients [38,39]. Moreover, LINC00538 was positively correlated with dendritic cells resting, in the meanwhile, negatively linked to dendritic cells activated in   our results. Previous studies demonstrated that the proportions of dendritic cells resting were positively correlated with the risk score of osteosarcomas [40] and LINC00538 was highly associated with worse outcomes of colon cancer patients [41]. ere are some disadvantages to our study as well. First, for starters, the study's findings need to be verified by bigger external partners. Second, more specific methods of interaction between m6A and lncRNAs, as well as how this regulation pattern contributes to the reshaping of TME, should be investigated further. ird, it should be examined if m6A-related lncRNAs participate in other biological processes associated with cancer. Lastly, while the associations between the risk score and histopathological characteristics or TME had a significant difference, the medical variations should be confirmed again because the variance was not as apparent. To validate the prediction models built in our current investigation, additional research studies should incorporate specimens from other databases as well as an increasing number of clinical specimens. More research is needed to completely understand the signaling pathways of m6Arelated lncRNAs in the carcinogenesis and development of osteosarcomas.
is study has some limitations: although the results of this study have some innovative and clinical significance, they have not been verified by more evidence, which makes this study incomplete, and a confirmatory study will be added in the future to make up for the shortcomings of this study.
In conclusion, our investigation showed that m6A-related lncRNAs remained tightly connected to the immunological microenvironment of osteosarcoma tumors, potentially influencing carcinogenesis and development. In addition, we looked for risk-related signaling of m6A-related lncRNAs in osteosarcomas and built and validated the risk prediction system. e findings of our current analysis will facilitate the assessment of the outcome and the development of immunotherapies for sufferers of osteosarcomas.
Data Availability e datasets used and analyzed during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.