Identified M6A-related LncRNAs in Patients with OS
The detailed flow diagram for the prognostic model is shown in Figure 1. Depending on the downloaded data file from the TCGA database, the matrix expression of 23 m6A genes and 14086 lncRNAs was screened. LncRNAs that were significantly related to greater than or equal to one of the 23 m6A genes (|Pearson R|༞0.4; P-value༜0.001) were defined as m6 A-related lncRNAs. Then, 535 m6A-related lncRNAs were obtained via Pearson correlation analysis and validated by the GEO database (GSE21257, GSE28424, and GSE19276). Finally, the m6A-lncRNA coexpression network was visualized using the network chart (Figure 2A). Given that lncRNAs exert a decisively favorable influence on tumorigenesis and development, univariate Cox regression was further performed to uncover the underlying relations between m6A-related lncRNAs and prognosis value (P value༜0.05). Results demonstrated that 43 m6A-related prognostic lncRNAs were significantly correlated with prognosis (Figure 2B). As well, the co-expression relationships between every m6A-related prognostic lncRNAs and immune checkpoints were also analyzed. Accordingly, SNHG12, LINC01357, SNHG7, GAS5, OLMALINC, ADAMTS9-AS2, SNHG6, AC079089.1, SNHG1, and SNHG4 were positively correlated with PD1 and CTLA-4, while FGD5-AS1, MSC-AS1, AC090559.1, and AC008074.2, were negatively correlated with them (Figure 2C and D). In summary, the results suggest that there are significant associations between immune checkpoints and m6A-related lncRNAs in OS.
Consensus Clustering of M6A-Related Prognostic LncRNAs in OS Patients
To uncover the underlying role of prognosis-associated m6A-related lncRNAs in the development of OS, consensus clustering was used to group OS patients, cluster1 and cluster2, based on the expression of 14 m6A-related prognostic lncRNAs. The matrix heat map for k = 2, also called CM plots, reveals the classification effect between the two clusters (Figure 3A). To find out the k for which the distribution reaches an approximate maximum, indicating maximum stability, the common clusters for k = 2 through 9 were reflected by the empirical cumulative distribution function (CDF) (Figure 3B). The score of the delta area plot illustrates the relative increase in cluster stability (Figure 3C). Together, a total of 85 patients with OS were divided into clusters 1 (n = 36) and 2 (n = 49).
To evaluate the different prognoses of patients between the two clusters, a Kaplan-Meier survival analysis was conducted. Results demonstrated that patients of cluster 2 had a longer survival time than cluster 1 (Figure 3D, P = 0.001). To figure out the possible associations between the m6a prognosis-related lncRNAs and clinical characteristics (metastasis, gender, and age), a heat map A heatmap was used to analyze their relationship (Figure 3E). Results demonstrated these clinical features were no significant difference in the two clusters. Meantime, the expression of most m6A-related lncRNAs such as AC110015.1, LINC01549, and AC010609.1 was highly expressed in cluster 1, while others such as PAXIP1−AS2 and MSC−AS1 were highly expressed in cluster 2.
Potential Biological Functions of the Two Clusters
To further clarify the potential biological process and pathway relating to the specific molecular differences between two subgroups, Gene Set Enrichment Analysis (GSEA) including Gene Ontology (GO), KEGG pathway, and tumor hallmark analysis was carried out from the TCGA database in different clusters (Table S1). In terms of GO cell component (CC), the gene sets were enriched in the basement membrane and collagen-containing extracellular matrix (ECM) (Figure 4A). The ECM structural constituent and receptor serine-threonine kinase binding are significantly enriched categories in GO molecular function (MF) (Figure 4B). The GO biological processes (BP) were mainly enriched in vitamin biosynthetic process, pigment granule localization, prostaglandin, and astrocyte development relative to cluster 2, while negative regulation of cell division and replacement ossification were enriched in cluster 1 (Figure 4C). Furthermore, KEGG analysis showed that the gene sets were mainly enriched in focal adhesion, mitogen-activated protein kinase (MAPK) signaling pathway, Hedgehog (Hh) signaling pathway, and so on, compared with cluster 1 (Figure 4D). Additionally, Hallmark analysis revealed that the gene sets were mainly enriched in apical surface, coagulation, apical junction, Hh signaling, epithelial-mesenchymal transition, and cholesterol homeostasis (Figure 4E). The outcome revealed that the cellular biological effects related to the m6A-related lncRNAs in OS.
M6A-related lncRNAs and TME in OS Patients
Herein, we further explore the immune infiltration ratio of different tumor-infiltrating immune cells in OS tissues. The CIBERSORT algorithm was used to analyze 22 different immune cell types in two clusters. Results revealed that infiltration fraction of CD8+ T cells, gamma delta T cells, resting NK cells, Monocytes, and M0 Macrophages showed obvious differences in two clusters (Figure 5A). Compared with cluster 2, cluster 1 has the lower levels of all immune cell types including CD8+ T cells, resting NK cells, and Monocytes, and higher levels of gamma delta T cells and M0 macrophages (Figure 5B). By the important roles of TME in tumorigenesis, the ESTIMATE evaluation is performed to examine stromal and immune scores for all OS samples. Here, there was no significant difference in immune scores between the two clusters. Furthermore, the ESTIMATE and stromal scores of cluster 2 were higher than that of cluster1 (Figure 5C). These results may give us some insights into the m6A-related lncRNA patterns that may remarkably influence specific immune cell types and TME, thus potentially suppressing or strengthening the response to immunotherapy.
Prognostic Analysis of Risk Model in M6A related lncRNAs
LASSO Cox regression analysis is extensively applied to heighten the forecast accuracy and interpretability of the statistical model. The application of this method is available to select the best prognostic markers to predict clinical results[14]. Consequently, 17 out of 43 m6A-related lncRNAs were discerned for the subsequent analysis (Figures 6A). To develop a signature for prognosis prediction of OS, patients (n=85) were divided into train set (n=44) and internal test set (n=41). Based on the median value of the prognostic risk grade, all sets were used to construct a high-low risk group to assess the prognostic risk of OS patients. The distribution of risk grade and the survival status of patients between the low-risk and high-risk groups are depicted in Figure 6B. Results revealed that the number of deaths is positively related to the risk score. Meanwhile, the relative expression standards of the 17 m6A-related lncRNAs for each patient are shown in Figure 6C. The expression of FGD5-AS1, SRP14-AS1, AC106771.1, ZFPM2-AS1, AC090559.1, AC009113.1, AC008074.2 were associated with the low-risk group, while BCAR3-AS1, AC110015.1, AC104461.1, AL121957.1, GAS5, MIR210HG, AC010609.1, AC025917.1, L161729.1, and SNHG4 were associated with the high-risk group. Kaplan-Meier survival analysis demonstrated that OS patients with the overall survival of higher risk scores were worse than lower risk scores (P༜0.001) (Figure 6D).
Validation of Risk Model in M6A related lncRNAs
To validate the predictive capability of this established model, the distribution of risk grade, survival status, and expression of the m6A-related lncRNAs were assessed in the test set (Figure 7A and B). Kaplan-Meier Survival analyses showed the overall survival of the low-risk group was longer than that of the high-risk group (P༜0.001) (Figure 7C). Results were consistent with the train set. Subsequently, we further analyzed the ROC curves in train and test sets (Figure 7D) and proved that candidate m6A-related lncRNAs had excellent predictive abilities for OS (AUC = 0.938 and 0.924, respectively). Finally, a heatmap for the correlated analysis was used among 17 m6A-related lncRNAs and the clinic characteristics (gender, age, and metastasis), clusters, and immune scores, respectively (Figure 7E and F). The expression of m6a prognosis-related lncRNAs in high-low risk groups was not different, while risk score was related to metastasis, immune scores, and clusters (P<0.05).
Independent prognostic analysis of Risk Model
Based on univariate and multivariate Cox regression analyses, the independent prognostic analysis of clinical features was performed to evaluate the risk model. We compared the risk scores and clinicopathological parameters including age, gender, and metastasis of OS in the train and test model, respectively (Figure 8). We found that risk score was an independent prognostic risk factor for OS (P<0.05). Furthermore, the subgroup analysis was applied to validate whether our model is useful to patients with different clinical parameters such as age, sex, and metastasis (Figure 9). Results revealed that our model could be used to the following different clinical characteristics: age (<=15 or >15), sex, and metastasis (P<0.05). These results have shown that the prognostic risk model of the 17 m6A-related lncRNAs for OS was comparatively dependable.
Evaluation of Risk Model with Immune cell substyles and immune checkpoint
To assess the relationship of immune cell substyles and immune checkpoints with risk scores, a scatter plot was depicted to validate this nexus. we found a positive correlation between risk scores and numbers of M0 macrophages, and negative correlations between risk scores and numbers of CD8 T cells and activated memory CD4 T cells (Figures 10A-C). The immune checkpoint expression of PD-1 and CTLA-4 were negatively correlated with the risk score (Figures 10D and E).