Identification and validation of an anoikis-related lncRNA signature to predict prognosis and immune landscape in osteosarcoma

Background: Anoikis is a specialized form of programmed apoptosis that occurs in two model epithelial cell lines and plays an important role in tumors. However, the prognostic value of anoikis-related lncRNA (ARLncs) in osteosarcoma (OS) has not been reported.

Methods: Based on GTEx and TARGET RNA sequencing data, we carried out a thorough bioinformatics analysis. The 27 anoikis-related genes were obtained from the Gene Set Enrichment Analysis (GSEA). Univariate Cox regression and least absolute shrinkage and selection operator (LASSO) analysis were successively used to screen for prognostic-related ARLncs. To create the prognostic signature of ARLncs, we performed multivariate Cox regression analysis. We calculated the risk score based on the risk coefficient, dividing OS patients into high- and low-risk subgroups. Additionally, the relationship between the OS immune microenvironment and risk prognostic models was investigated using function enrichment, including Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), single-sample gene set enrichment analysis (ssGSEA), and GSEA analysis. Finally, the potential effective drugs in OS were found by immune checkpoint and drug sensitivity screening.

Results: A prognostic signature consisting of four ARLncs (AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1) was constructed. To assess the regulation patterns of anoikis-related lncRNA genes, we created a risk score model. According to a survival analysis, high-risk patients have a poor prognosis as they progress. By using immune functional analysis, the lower-risk group demonstrated the opposite effects compared with the higher-risk group. GO and KEGG analysis showed that the ARLncs pathways and immune-related pathways were enriched. Immune checkpoints and drug sensitivity analysis might be used to determine the better effects of the higher group.

Conclusion: We identified a novel prognostic model based on a four-ARLncs signature that might serve as potential prognostic indicators that can be used to predict the prognosis of OS patients, and immunotherapy and drugs that may contribute to improving the overall survival of OS patients and advance our understanding of OS.

Introduction

Osteosarcoma (OS) is the most common primary malignant bone tumor and is common in children and adolescents (1). It has an annual prevalence of three to four cases per million people, but it is an uncommon condition (2). Neoadjuvant chemotherapy is linked to a better 5-year overall survival in patients with limb osteosarcoma (3). Osteosarcoma generally develops in weight-bearing long bones, with the distal femur and proximal tibia being the most frequent sites (4). Compared with surgery and chemotherapy, a standard therapy for osteosarcoma has been developed; however, therapeutic advancement is still limited (5, 6). Osteosarcoma’s poor prognosis was strongly correlated with both its high recurrence rate and distant lung metastasis (7). Numerous studies built on clinical knowledge are looking into osteosarcoma therapeutic approaches. However, finding precise and efficient therapy targets for OS prognosis is challenging (8). Therefore, in order to develop new therapeutic strategies, it is crucial to investigate the regulatory mechanisms of genes related to osteosarcoma. lncRNAs, also known as transcripts, are ncRNAs that are longer than 200 nucleotides and do not code for proteins (9–11). Long non-coding RNAs (lncRNAs) participate in multiple biological processes and regulate physiological functions (12). Anoikis is a specialized form of programmed apoptosis that occurs in two model epithelial cell lines when cell adherence to the ECM is lost or cell–matrix interactions are disrupted (13, 14). Anoikis resistance is a property that cancer cells can acquire, allowing them to grow and spread from their primary sites to distant organs (15). By preventing detached cells from adhering to other substrates for abnormal proliferation, apoptosis plays a crucial role in protecting the organism (16). Numerous studies have demonstrated a connection between anoikis-related genes (ARGs) and the tumor metastasis cascade and cancer development (17)—for instance, human osteosarcoma cells become anoikis resistant by upregulating VEGF-A expression, which may promote angiogenesis and lung metastasis in OS (18). Numerous studies have revealed that lncRNA dysregulation is a significant factor in the modulation of anoikis resistance in a variety of cancers (19). lncRNAs have the ability to positively or negatively influence anoikis resistance (20). There are not many papers that analyze prognostic markers based on ARLncs in OS, even though many cancers have been shown to be related to anoikis (21–23). Thus, it is significant to focus on the relationship between integrated ARGs and OS clinical outcomes to investigate the therapeutics of OS. Using prognostic models, genes associated with necroptosis and apoptosis have been found to predict OS prognosis (24, 25). In this study, we developed a novel OS prognostic model to investigate anoikis-related lncRNA (ARLncs), which can be used to develop risk models and immunological features. We identified a robust ARLncs-based signature and exploited its clinical implications in OS patients. It is helpful for us to explore the prediction of how the tumor immunity microenvironment is regulated by the ARLncs signature, which has promising application possibilities.

Methods

Data acquisition

The gene expression and clinical information of 88 osteosarcoma patients were taken from the TARGET database (https://portal.gdc.cancer.gov/). From the GTEx database (https://gtexportal.org/), the gene expression information for 803 samples of normal muscular tissues was gathered. We excluded the samples that lacked sufficient information and divided the remaining samples into training and test groups of 42 each using a 1/1 ratio analysis. The expression data of gene for each sample were transformed into log2 (x + 1) to remove data differences caused by different platforms between the TCGA and GTEx databases. Batch-to-batch variation was then eliminated by combining two datasets into a single dataset using the combat function from the “sva” R package.

Differentially expressed anoikis-related lncRNAs

With |log2FC| >1 and false discovery rate (FDR) <0.05, we use the R package “limma” to identify the lncRNA genes that are differentially expressed between 88 OS tissues and 803 normal tissues. The 27 anokis-related genes were collected from the gene set enrichment analysis (GSEA). According to Pearson’s correlation, the anoikis-related lncRNAs were found with Cor >0.4 and P <0.001 (26).

Construction of anoikis-related lncRNA prognostic signatures

The “survival” package in R was used to identify lncRNAs associated with patient prognosis among the anoikis-related DE-lncRNAs by using uni-Cox regression with a P-value <0.05. To identify anoikis-related lncRNAs, we utilized 10-fold cross-validation and LASSO regression. Additionally, the anoikis-related lncRNA signature was chosen as the modeling gene after performing a multi-Cox regression analysis using the “glmnet” package to identify the optimum penalty value. The following formula calculates the risk score of OS patients:

$$\text{Risk score }= \sum _{\text{i}=1}^{\text{n}}({\text{lncRNAexp}}_{\text{i}} \times {\text{coef}}_{\text{i}})$$

By computing the median risk score according to the algorithm, the OS patients were divided into low- and high-risk groups. The KM survival analysis of anoikis-related lncRNA signatures was built in the whole and in subgroups by the “survival” package in R software. Based on the timeROC R packages, the time-related and clinical-related ROC curves were built to investigate the predictive capacity of this gene signature. To determine if risk score can be seen as a factor influencing the prognosis of OS patients, the survival package in R was employed (27).

Construction of nomogram, calibration, clinical forest, C-index, PCA, and t-SNE

Combined with risk scores, the nomogram was set up with the rms R package based on gender, age, tumor site, and metastasis. The C-index was established to evaluate the prediction bias of the nomogram. According to related R packages, we combined clinical data and risk scores to build a forest tree by uni- and multi-Cox regression analysis. To assess the efficacy of the risk model in OS patients in the low- and high-risk groups, principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were constructed, respectively, on the risk prognosis in the entire, training, and test samples.

Functional analysis, GSEA analysis, and immunological microenvironment analysis

In order to examine the underlying biological process, we determined the differential genes between the low-risk group and the high-risk group for volcano plot, heat map, and later Gene Ontology (GO) analysis. Additionally, GSEA was used to determine the KEGG pathways that were differentially expressed in the two groups. The KEGG gene set (c2.cp.kegg.v2022.1.Hs.symbols.gmt) can be obtained from the website (https://www.gsea‐msigdb.org/). The website served as the source for the KEGG gene set. P-value <0.05 was chosen as the cutoff for biological processes and pathways that were significantly enriched. Using R software’s limma, reshape2, and ggpubr packages, we examined the expression differences of immune cells and immune function between the low- and high-risk groups of the entire sample (28).

Analysis of immune checkpoints and drug sensitivity

The “pRRophetic” package of R is commonly used to perform chemotherapeutic prediction (29, 30). The R package “pRRophetic” is utilized to investigate the difference in drug sensitivity among OS patients. Drug sensitivity analysis, using the standard of P <0.001, was used to identify whether the sensitivity of drugs was different in the risk prognosis model including all samples. We evaluate the expression of immune checkpoint-related genes in the low- and high-risk groups of OS patients to examine the impact of immunotherapy with the help of “reshap2”, “limma”, “ggplot2”, and “ggpubr” packages.

Results

Identification of a risk prognostic model 8in osteosarcoma

We conducted differential analysis to get 1,214 lncRNA-related differentially expressed genes between 88 OS tissues and 803 normal tissues with |log2FC| >1 and FDR <0.05. Then, we carried out Pearson’s correlation analysis to confirm 184 DEAKLncs with 27 anoikis-related genes from the GSEA database. Anokis-related genes and lncRNA were correlated and integrated into a network as shown in Figures 1A, B. First, 84 OS samples were randomly divided into two groups in a ratio of 1:1. Then, we constructed univariate Cox regression. The samples were divided into two groups based on the risk model, with more information about each group provided in Table 1. Figure 1C displays a heat map between OS tumor tissues and normal tissues and a nine-gene forest tree by univariate Cox regression (Figure 1D). The optimal value of ARLncs participating in the model construction by LASSO regression analysis was 6 (Figures 1E, F). Through multi-Cox regression, the forest tree in Figure 1G was revealed. Lastly, four target genes (AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1) were identified through uni-Cox regression, LASSO regression, and multi-Cox regression analyses to construct the prognostic model of both risk groups in OS patients, with three lncRNAs (MEF2C-AS1, SNHG6, and TBX2-AS1) considered as risk factors and one lncRNA (AC079612.1) as a protective factor for osteosarcoma patients. The coefficients for the four target genes can be found in Table 2. Through the formula for calculating the risk score, we can get the risk score for all samples by calculating them as follows: risk score = (-0.422 * AC079612.1exp) + (0.963 * MEF2C-AS1exp) + (0.619 * SNHG6exp) + (0.751 * TBX2-AS1exp).

FIGURE 1

Figure 1 Developing a four anoikis-related long noncoding RNA (lncRNA) signature for osteosarcoma patients. (A) Sankey diagram highlighting the relationship between the anoikis genes and the anoikis-related lncRNAs. (B) Co-expression network showing the relationship between the expression of lncRNAs and that of Anoikis-associated genes. (C) Forest plot highlighting the univariate Cox regression analysis results. (D) Heatmap showing differential expression of Anoikis-related lncRNAs in tumor tissues and normal tissues. Red indicates high expression levels, while blue indicates low expression levels. (E) A general cross- validation curve of the paired likelihood deviance. (F) Elucidation for the LASSO coefficient profiles of the prognostic lncRNAs. (G) Forest plot highlighting the multiple Cox regression analysis results.

TABLE 1

Table 1 Clinical features of all osteosarcomas categorized into the training and validation cohorts.

TABLE 2

Table 2 Optimal prognostic risk signature of four lncRNA by multivariate Cox regression analysis.

We compared the expression of four target anoikis-related lncRNAs in the low- and high-risk groups to evaluate the prognostic capability of the risk model, the distribution of risk score, and the survival status of the four ARLncs between risk subgroups in the whole group, training group, and test group (Figures 2A–C). In a Kaplan–Meier survival study of three groups, the prognosis of the low-risk group was statistically significantly better than that of the high-risk group. Figure 3A displays the relationship between the four strong prognostic lncRNAs and the 27 anoikis genes throughout the whole TARGET OS cohort.

FIGURE 2

Figure 2 Prognosis of the risk model in the entire, training, and testing sets. (A–C) Survival time, survival status, Kaplan–Meier survival curves of patients with OS, and heat map in the expression of four lncRNAs in the entire, training, and testing sets, respectively.

FIGURE 3

Figure 3 Combination of four DEARLncs signature and clinical features, analysis of PCA and t-SNE. (A) Relationship between four target genes and Anoikis-related lncRNAs in the entire cohort. (B) Nomogram based on age, gender, tumor site, metastasis, and risk in the TARGET database. (C) Calibration curves illustrated the consistency between predicted and observed 1-, 3-, and 5-year survival rates in OS patients depending on the prognostic nomogram. (D, E) The ROC curves for 1, 3, 5 years and clinical survival in the TARGET cohort. (F) The C-index was used to determining the prognostic accuracy of risk scores and addition clinical factors in entire cohort. (G–I) The PCA and t-SNE in entire, training, testing sets, respectively.

The prognosis of OS patients was predicted using a nomogram that contained clinical–pathological factors and risk scores at the 1-, 3-, and 5-year intervals (Figure 3B). According to the calibration curves, the actual overall survival rates and the anticipated survival rates at 1, 3, and 5 years were in good agreement (Figure 3C). The AUCs of the 1-, 3-, and 5-year survival rates for the ARLncs signature were 0.872, 0.774, and 0.807, respectively (Figure 3D). When combined with clinical data, the AUCs of risk, gender, age, metastasis, tumor site, and the AUCs of risk were the best in predicting overall survival in OS patients; metastasis ranked second (Figure 3E). The C-index was used to determining the prognostic accuracy of risk scores and addition clinical factors in entire cohort (Figure 3F). The expression levels of the OS prognostic ARLncs were used in the risk model construction of the entire, training, and test cohorts, respectively, and were shown to significantly identify patients in the low- and high-risk groups, which can demonstrate the accuracy of the model. This was displayed by PCA and t-SNE (Figures 3G–I). Above all, the prognostic model can be well predicted in osteosarcoma patients.

Independent prognostic signature

The overall survival-related variable between risk score and clinical factors was analyzed using the univariate Cox regression method (Figure 4A). The outcome demonstrated a statistically significant difference between metastatic status and risk score. The risk score and metastatic status were found to be independent risk factors affecting the prognosis of OS when we utilized multivariate Cox regression to analyze the risk score and clinical factors at the same time (Figure 4B). In order to compare the survival differences of the low- and high-risk groups in various clinical subgroups, we divided the OS patients into groups based on gender, age (≤14 and >14), and metastasis. The results showed that OS patients had a better clinical prognosis in the low-risk group than in the high-risk group among different genders, ages, and metastatic statuses (Figures 4C–E). Therefore, the risk model constructed by the four anoikis-related lncRNAs in the low-risk group would affect the survival prognosis of OS patients, which, combined with metastasis, would be a novel prognosis signature. The expression levels of anoikis-related lncRNAs (AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1) showed differences in risk subgroups (Figures 5A–D) and different tissue types (Figures 5E–I). In the risk subgroups, the expression level of AC079612.1 was lower in the high-risk group than in the low-risk group, acting as a protective factor (Figure 5A), but other genes (MEF2C-AS1, SNHG6, and TBX2-AS1) were higher in the high-risk group (Figures 5B–D). In contrast, the higher differential expression level of AC079612.1 in tumor tissues was obvious compared with normal tissues (Figure 5E).

FIGURE 4

Figure 4 Prognosis of the anoikis-related lncRNA signature in osteosarcoma patients. (A) Univariate COX regression analysis. (B) Multivariate COX regression analysis. (C) Survival curve of gender. (D) Survival curve of age. (E) Survival curve of metastasis.

FIGURE 5

Figure 5 Differential expression of target genes in different osteosarcoma risk and tumor and normal tissues for the expression of AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1 in various risk groups (A–D) and different tissue types (E–I). "ns" P > 0.05, *** P < 0.001.

Functional enrichment analysis

Differentially analyzing the genes in both risk groups helped identify the potential molecular mechanisms of the four-anoikis-related lncRNA signatures. A total of 833 DEGs were found with |logFC| >1 and a P-value <0.05. In total, 322 genes were upregulated, and 501 genes were downregulated among the DEGs. Figures 6A, B, respectively, display the volcano plots and heat map of anoikis-related lncRNAs with differential expression in the low- and high-risk groups. A GO analysis that we conducted revealed that the DEGs were primarily enriched in immune activation and immune response pathways. The GO bubble displayed a positive regulation of leukocyte activation in BP, exterior side of plasma membrane in CC, and antigen binding in MF (Figure 6C). Figure 6D shows the different function ratio of GO terms in GO circle. The GO chord showed that the positive regulation of lymphocyte activation and the positive regulation of cell activation were mainly enriched (Figure 6E). In the GO analysis, the network of biological processes is displayed in Figure 6F. According to the KEGG analysis, the top 22 KEGG terms in Figures 7A, B were primarily enriched in the cytokine–cytokine receptor interaction and PI3K-Akt signaling pathway. Overall, GO and KEGG analysis showed a strong link between anoikis and immunological status. To determine the potential process and pathway between the low- and high-risk groups, GSEA was carried out for the KEGG pathway enrichment analysis (Figures 7C, D). Most immune-related pathways, including antigen processing and presentation, cytokine–cytokine receptor interaction, hematopoietic cell lineage, Nod-like receptor signaling pathway, and primary immunodeficiency, were proven to exist in the low-risk group. On the other hand, several pathways were involved in OS development, including insulin signaling pathway, long-term potentiation, ribosome, and terpenoid backbone biosynthesis, and were mainly enriched in the high-risk group. Overall, the results showed that these tumor and immune-related pathways were regulated by the novel anoikis-related lncRNA signature. The microenvironment and immune infiltration of tumor in both risk groups

FIGURE 6

Figure 6 Differential analysis and GO functional enrichment analysis of the TARGET database in the high- and low-risk groups. (A) Volcano plot of differentially express(DEGs): red for high expression and blue for low expression. (B) Heat map of DEGs, with high expression in red and low expression in green. (C) Gene Ontology (GO) analysis for the DEGs between the high- and low-risk groups. (D) Circle plots for the analysis. (E) GO enrichment analysis of DEARLncs. (F) Network of biological processes in GO analysis.

FIGURE 7

Figure 7 Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Set Enrichment Analysis (GSEA) of different risk subgroups. (A) Bar plot of KEGG. (B) Bubble of KEGG. (C) Multiple pathways and functions were found to be enriched in low-risk anoikis-related lncRNA. (D) GSEA analysis in high-risk groups.

We investigated the infiltration of 22 immune cell types in OS patients using the CIBERSORT algorithm. Most immune cell filtration were naive CD4 T cells, macrophages M0, M1, and M2, and activated dendritic cells (Figure 8A). According to risk groups, the expression of naïve CD4 T cells was higher in the high-risk group than in the low-risk group (P = 0.004); other immune cells did not have statistical significance in both risk groups (Figure 8B). The levels of immune cell infiltration were evaluated using the single sample gene set enrichment analysis (ssGSEA) based on the TARGET dataset (31). The box plot revealed that immune cell infiltration had lower scores in high-risk groups than in low-risk groups. However, the aDC infiltration had a higher score in the high-risk subgroup (Figure 8C). The most functional pathways including APC_co_inhibition, Check-point, Cytolytic_activity, T_Cell_co-inhibition, and T_Cell_co-stimulation were mainly enriched in the low-risk group (Figure 8D). We further analyzed the relationship between immune cell infiltration and risk score to investigate the function of anoikis-related lncRNAs in the OS immune microenvironment. We evaluated the ESTIMATEScore, ImmuneScore, and StromalScore, which were higher in the low-risk group than in the high-risk group. Meanwhile, the scores showed a negative correlation with riskScore (Figures 8E–G). In contrast to the high-risk group, TumorPurity was lower in the low-risk group. Similar to the box plot, the correlation was positive with riskScore (Figure 8H). When combined with clinical data, the heat map showed that the expression of three anoikis-related lncRNAs (MEF2C-AS1, SNHG6, and TBX2-AS1) was higher in the high-risk group compared with the low-risk group, while AC079612.1 was expressed at a lower level in the high-risk group. At the same time, statistical significance (P < 0.001) revealed that the high-risk group had a strong correlation with osteosarcoma metastasis (Figure 9A). The heat map shows the infiltration of immune cells in different risk groups combined with the immune microenvironment—for example, the expression of Type_II_IFN_Response was higher in the high-risk group, and there was a higher score in Stromalscore, ImmuneScore, ESTIMATEScore, and TumorPurity (Figure 9B). We know that the immune-related functions were significantly downregulated in the high-risk group (Figure 9C). Regarding the correlation between tumor microenvironment infiltration and riskScores, ImmuneScore, Stromalscore, ESTIMATEscore, and TumorPurity in the risk groups, we could see that riskScore was related with B cell memory, activated dendritic cells, Stromalscore, ESTIMATEscore, and TumorPurity (Figure 9D). The pathways in GSVA were regulated by the AC079612.1 gene, such as the pathways including kegg_circadian_rhythm_mammal and kegg_ascorebate_and_aldarate_metabolism, which were enriched in the Up group (Figure 9E).

FIGURE 8

Figure 8 Immune signature in the low- and high-risk groups. (A) Proportion of 22 immune cell types in osteosarcoma risk groups. (B) Expression of immune cells between two groups. (C) The differential expression of cell infiltration between risk groups is based on the ssGSEA scores. (D) Immune functional difference between risk groups. (E–H) Differential expression of tumor microenvironment scores (immune scores, stromal scores, ESTIMATE scores, and tumor purity) between risk groups. * P < 0.05, ** P < 0.01, *** P < 0.001.

FIGURE 9

Figure 9 Distribution and visualization of immune status. (A) Heat map and the clinical characters of the two groups (***P < 0.001). (B) Heat map for immune status based on ESTIMATE and ssGSEA among two risk subgroups. (C) Heat map visualizing the infiltrated immune cells in the low- and high-risk groups. (D) Correlation between immune cells and riskScores, ImmuneScore, StromalScore, ESTIMATEScore, and TumorPurity. *P < 0.05, **P < 0.01, ***P < 0.001. (E) Count of pathway by AC079612.1 gene in GSVA.

Correlation between immune cells and target genes, drug sensitivity, and immune checkpoint molecule expression analysis

Under a statistical significance, Dendritic cells resting and Macrophages M1 had a positive correlation with AC079612.1, and T cells regulatory (Tregs) was shown to be negatively correlated with AC079612.1. Resting dendritic cells were positively correlated with MEF2C−AS1 but had a negative correlation with Tregs (Figure 10). Based on the above-mentioned data, AC079612.1 gene might be a good biomarker to predict the therapy and diagnose osteosarcoma. Therefore, we explored the correlation between AC079612.1 and 27 anoikis-related genes by conducting Pearson’s correlation analysis. The results showed that AC079612.1 expression was positively correlated with the expression of PDK4 and SRC genes (|R| > 0.4, P < 0.05) (Figures 11A–Z2). The drug sensitivity analysis revealed that AZD7762, bortezomib, etoposide, gemcitabine, KIN001-135, MP470, T0901317, temsirolimus, and thapsigargin (P < 0.001) all showed statistically significant sensitivity in both risk groups. AZD7762, bortezomib, and etoposide had more drug sensitivity in both risk groups, while patients in the high-risk group were more sensitive to gemcitabine, KIN001-135, MP470, T0901317, temsirolimus, and thapsigargin than those in the low-risk group (Figures 12A–I). The different drugs were sensitive in the treatment of OS patients by drug sensitivity analysis. Given the importance of checkpoint-based immunotherapy, we noted that the expression of TNFRSF9, BTLA, CD200R1, TIGIT, HAVCR2, CD274, PDCD1LG2, TNFSF15, NRP1, CD44, CD27, LAG3, CD276, LGALS9, LAIR1, KIR3DL1, and CD40LG was higher in the low-risk group than in the high-risk group (Figure 12J), showing that OS patients could better respond to immunotherapy.

FIGURE 10

Figure 10 Correlation between immune cells (Dendritic.cells, T.cells.regulatory, Macrophages.M1, T.cells.regulatory, and Dendritic.cells) and target genes (AC079612.1 and MEF2C-AS1).

FIGURE 11

Figure 11 Pearson correlation analysis of AC079612.1 with anoikis-related genes (A–Z2: CRYBA1, ITGB1, MAP3K7, NOTCH1, CEACAM6, IKBKG, CAV1, MCL1, BRMS1, SNAI2, ITGA5, MTOR, DAPK2, PIK3CA, PTRH2, PDK4, TLE1, TSC2, E2F1, CHEK2, CEACAM5, PTK2, AKT1, NTRK2, STK11, SRC, and BCL2).

FIGURE 12

Figure 12 Drug sensitivity analysis and immune checkpoint molecules expression analysis. (A–I) Sensitivity performance of 9 drugs in the high-risk and low-risk groups. (J) Differential expression of immune checkpoint genes between low- and high-risk groups. * P < 0.05, ** P < 0.01.

Discussion

The majority of occurrences of osteosarcoma, the most common malignant bone tumor, are seen in children and young people between the ages of 10 and 30 (32). Despite the fact that standard treatments, such as a series of comprehensive treatments including surgery and chemotherapy, increased the survival rates, certain patients’ odds of surviving with OS were not good (33). Drug resistance during the therapy process and OS metastases caused the treatment effects to be unsatisfactory (34). Therefore, it is imperative to find reliable prognostic indicators for treating OS patients.

Anoikis, which prevents shed cells from sticking to a new substrate in the wrong location and restricting the organism’s growth, is an essential defense mechanism (35). According to certain studies, the intrinsic and extrinsic mechanisms must both be present for anoikis apoptosis to take place (36). While mitochondria are important to regulate apoptosis, a number of intracellular signals, such as DNA damage, can also cause anoikis (37). This disorder in anoikis execution, which promotes tumor cell invasion and migration, the organ’s growth through metastasis, and the establishment of drug resistance, may be present in cancer cells (38–40). There were increasing studies that showed anoikis resistance and pulmonary metastasis in the OS microenvironment by altering related signal pathways (41).

Additionally, a growing amount of data suggests that lncRNA plays critical roles in controlling the development of various cancers (42). The prognosis of tumor patients was predicted using a variety of novel lncRNA models based on anoikis (19, 43). Currently, there is no research to report on the study of both anoikis-related lncRNA and the clinical characteristics of OS patients. Hence, finding the anoikis-related lncRNAs as the prognosis markers of osteosarcoma is essential. In this study, we conducted a comprehensive analysis of the target lncRNAs and the related immune microenvironment to explore the potential molecular mechanisms of OS progression.

In this study, the expression levels of the ARGs in OS and normal tissues were investigated to assess the prognostic efficacy of anoikis-related lncRNAs in OS patients. Then, the uni-Cox regression analysis was used to identify nine prognostic ARLncs. Four important ARLncs, including AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1, were screened using the LASSO Cox and Multiple Cox regression analysis for building the optical model, and a four-ARLncs signature for osteosarcoma was effectively constructed. Of the four anoikis-related lncRNAs, we discovered that AC079612.1 was a protective factor with a lower expression in the high-risk subgroup of osteosarcoma. Patients at high risk have a significantly lower rate of survival than patients at low risk. The validation cohort’s prognosis also demonstrated the model was well. Additionally, a great nomogram for predicting survival rates was developed using the risk scores and additional clinical factors (such as gender, age, metastasis, and tumor site). In conclusion, these results demonstrated the four-ARLncs signature’s applicability to other cohorts and supported its strong predictive significance in our study’s osteosarcoma patients. The tumor microenvironment, which played a significant role in the treatment of OS, was identified as one of the most critical factors affecting immunotherapy (44, 45). Currently, immune-related therapy has become an increasingly important strategy for OS patients (46, 47). We explored functional enrichment analysis using GO, KEGG, and GSEA to investigate the anoikis-related lncRNAs in OS patients. The pathways of antigen processing and presentation, cytokine–cytokine receptor interaction, hematopoietic cell lineage, and Nod-like receptor signaling pathway were found in related immune-mechanisms of OS. Next, the relationship between the anoikis-related lncRNAs signature and the immune status of osteosarcoma was then analyzed. The tumor microenvironment in OS was explored using the ESTIMATE, ssGSEA, and CIBERSORT algorithms. Finally, we identified the correlation between risk scores and immune cells and the association of target ARLncs with anoikis-related genes. Of course, drug sensitivity and immune checkpoint studies were conducted to explore the effectiveness of drug therapy and immunotherapy in both risk groups.

Anoikis-related lncRNAs are closely related to the mechanism of tumorigenesis, but the relationship between osteosarcoma and ARLncs is not known. Therefore, it is of great potential value to study the mechanism and treatment of osteosarcoma from the point of view of ARLncs. Through comprehensive analysis, four ARLncs genes and related immune pathways were identified, which is of great help in understanding the mechanism of osteosarcoma at the genetic level. The drug sensitivity analysis provides a new theoretical support for clinicians treating osteosarcoma with drugs.

However, there are some limitations in this study. First, due to the gene expression derived from the public TARGET database, the samples of osteosarcoma are small. Future work should focus on developing larger samples to improve the accuracy of the results. Second, there were only 84 clinical data added to the database; if possible, the clinical data can be collected in hospitals to increase the abundance of clinical data, and the data on tumor stage in OS was lacking. Third, the other testing dataset is needed to assess the value of the risk score model. Fourth, there is absence of a corresponding experiment with target genes to verify the expression and effectiveness of genes in vivo in clinical trials. Finally, the underlying mechanism of four ARLncs (AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1) regulating the prognosis of OS patients should be explored further. These drawbacks will be included in our follow-up plans and addressed. For the small sample size of osteosarcoma, we will collect clinical samples for further research in the future. At the same time, if conditions permit, we will also collect osteosarcoma samples for experimental verification, thereby increasing the reliability of the conclusion.

Conclusion

The risk model would be constructed using four ARLncs (AC079612.1, MEF2C-AS1, SNHG6, and TBX2-AS1). We found that the four ARLncs can influence OS’s immunologic microenvironment. Nine drugs and immunotherapy can be identified as potential therapeutical schemes for OS. These results of the study are helpful to increase the therapeutic effectiveness and improve the general survival rate of individuals with osteosarcoma.

Data availability statement

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 below: TCGA, TARGET-OS, phs000218.

Author contributions

J-SZ constructed the scheme and collected the data of the study. J-SZ drafted the manuscript, and R-SP revised the manuscript. Funding acquisition was from X-BT. All authors contributed to the article and approved the submitted version.

Funding

The study was funded by the Science and Technology Foundation of Guizhou Province (grant number GZKJ-2021-072).

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Rojas GA, Hubbard AK, Diessner BJ, Ribeiro KB, Spector LG. International trends in incidence of osteosarcoma (1988-2012). Int J Cancer (2021) 149(5):1044–53. doi: 10.1002/ijc.33673

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Smeland S, Bielack SS, Whelan J, Bernstein M, Hogendoorn P, Krailo MD, et al. Survival and prognosis with osteosarcoma: outcomes in more than 2000 patients in the EURAMOS-1 (European and American osteosarcoma study) cohort. Eur J Cancer (2019) 109:36–50. doi: 10.1016/j.ejca.2018.11.027

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Papakonstantinou E, Stamatopoulos A, Athanasiadis DI, Kenanidis E, Potoupnis M, Haidich AB, et al. Limb-salvage surgery offers better five-year survival rate than amputation in patients with limb osteosarcoma treated with neoadjuvant chemotherapy. a systematic review and meta-analysis. J Bone Oncol (2020) 25:100319. doi: 10.1016/j.jbo.2020.100319

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Shoaib Z, Fan TM, Irudayaraj JMK. Osteosarcoma mechanobiology and therapeutic targets. Br J Pharmacol (2022) 179(2):201–17. doi: 10.1111/bph.15713

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gill J, Gorlick R. Advancing therapy for osteosarcoma. Nat Rev Clin Oncol (2021) 18(10):609–24. doi: 10.1038/s41571-021-00519-8

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Isakoff MS, Bielack SS, Meltzer P, Gorlick R. Osteosarcoma: Current treatment and a collaborative pathway to success. J Clin Oncol (2015) 33(27):3029–35. doi: 10.1200/jco.2014.59.4895

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Huang X, Zhao J, Bai J, Shen H, Zhang B, Deng L, et al. Risk and clinicopathological features of osteosarcoma metastasis to the lung: A population-based study. J Bone Oncol (2019) 16:100230. doi: 10.1016/j.jbo.2019.100230

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Liu Y, Huang N, Liao S, Rothzerg E, Yao F, Li Y, et al. Current research progress in targeted anti-angiogenesis therapy for osteosarcoma. Cell Prolif (2021) 54(9):e13102. doi: 10.1111/cpr.13102

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Panni S, Lovering RC, Porras P, Orchard S. Non-coding RNA regulatory networks. Biochim Biophys Acta Gene Regul Mech (2020) 1863(6):194417. doi: 10.1016/j.bbagrm.2019.194417

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Peng WX, Koirala P, Mo YY. LncRNA-mediated regulation of cell signaling in cancer. Oncogene (2017) 36(41):5661–67. doi: 10.1038/onc.2017.184

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sanchez Calle A, Kawamura Y, Yamamoto Y, Takeshita F, Ochiya T. Emerging roles of long non-coding RNA in cancer. Cancer Sci (2018) 109(7):2093–100. doi: 10.1111/cas.13642

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang P, Wu W, Chen Q, Chen M. Non-coding RNAs and their integrated networks. J Integr Bioinf (2019) 16(3):20190027. doi: 10.1515/jib-2019-0027

CrossRef Full Text | Google Scholar

13. Frisch SM, Francis H. Disruption of epithelial cell-matrix interactions induces apoptosis. J Cell Biol (1994) 124(4):619–26. doi: 10.1083/jcb.124.4.619

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Wang J, Luo Z, Lin L, Sui X, Yu L, Xu C, et al. Anoikis-associated lung cancer metastasis: Mechanisms and therapies. Cancers (Basel) (2022) 14(19):4791. doi: 10.3390/cancers14194791

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Kim YN, Koo KH, Sung JY, Yun UJ, Kim H. Anoikis resistance: an essential prerequisite for tumor metastasis. Int J Cell Biol (2012) 2012:306879. doi: 10.1155/2012/306879

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Chen S, Gu J, Zhang Q, Hu Y, Ge Y. Development of biomarker signatures associated with anoikis to predict prognosis in endometrial carcinoma patients. J Oncol (2021) 2021:3375297. doi: 10.1155/2021/3375297

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Ye G, Yang Q, Lei X, Zhu X, Li F, He J, et al. Nuclear MYH9-induced CTNNB1 transcription, targeted by staurosporin, promotes gastric cancer cell anoikis resistance and metastasis. Theranostics (2020) 10(17):7545–60. doi: 10.7150/thno.46001

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Gao Z, Zhao GS, Lv Y, Peng D, Tang X, Song H, et al. Anoikis−resistant human osteosarcoma cells display significant angiogenesis by activating the src kinase−mediated MAPK pathway. Oncol Rep (2019) 41(1):235–45. doi: 10.3892/or.2018.6827

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lee HY, Son SW, Moeng S, Choi SY, Park JK. The role of noncoding RNAs in the regulation of anoikis and anchorage-independent growth in cancer. Int J Mol Sci (2021) 22(2):627. doi: 10.3390/ijms22020627

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Shi T, Zhang C, Xia S. The potential roles and mechanisms of non-coding RNAs in cancer anoikis resistance. Mol Cell Biochem (2022) 477(5):1371–80. doi: 10.1007/s11010-022-04384-6

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Jin L, Chun J, Pan C, Kumar A, Zhang G, Ha Y, et al. The PLAG1-GDH1 axis promotes anoikis resistance and tumor metastasis through CamKK2-AMPK signaling in LKB1-deficient lung cancer. Mol Cell (2018) 69(1):87–99.e7. doi: 10.1016/j.molcel.2017.11.025

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Simpson CD, Anyiwe K, Schimmer AD. Anoikis resistance and tumor metastasis. Cancer Lett (2008) 272(2):177–85. doi: 10.1016/j.canlet.2008.05.029

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wang YN, Zeng ZL, Lu J, Wang Y, Liu ZX, He MM, et al. CPT1A-mediated fatty acid oxidation promotes colorectal cancer cell metastasis by inhibiting anoikis. Oncogene (2018) 37(46):6025–40. doi: 10.1038/s41388-018-0384-z

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hong J, Li Q, Wang X, Li J, Ding W, Hu H, et al. Development and validation of apoptosis-related signature and molecular subtype to improve prognosis prediction in osteosarcoma patients. J Clin Lab Anal (2022) 36(7):e24501. doi: 10.1002/jcla.24501

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Wang G, Zhang X, Feng W, Wang J. Prediction of prognosis and immunotherapy of osteosarcoma based on necroptosis-related lncRNAs. Front Genet (2022) 13:917935. doi: 10.3389/fgene.2022.917935

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Fan Q, Liu B. Identification of a RNA-seq based 8-long non-coding RNA signature predicting survival in esophageal cancer. Med Sci Monit (2016) 22:5163–72. doi: 10.12659/msm.902615

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jiang F, Miao XL, Zhang XT, Yan F, Mao Y, Wu CY, et al. A hypoxia gene-based signature to predict the survival and affect the tumor immune microenvironment of osteosarcoma in children. J Immunol Res (2021) 2021:5523832. doi: 10.1155/2021/5523832

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhang Y, He R, Lei X, Mao L, Jiang P, Ni C, et al. A novel pyroptosis-related signature for predicting prognosis and indicating immune microenvironment features in osteosarcoma. Front Genet (2021) 12:780780. doi: 10.3389/fgene.2021.780780

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature (2009) 462(7269):108–12. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Meltzer PS, Helman LJ. New horizons in the treatment of osteosarcoma. N Engl J Med (2021) 385(22):2066–76. doi: 10.1056/NEJMra2103423

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Durfee RA, Mohammed M, Luu HH. Review of osteosarcoma and current management. Rheumatol Ther (2016) 3(2):221–43. doi: 10.1007/s40744-016-0046-y

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Hameed M, Dorfman H. Primary malignant bone tumors–recent developments. Semin Diagn Pathol (2011) 28(1):86–101. doi: 10.1053/j.semdp.2011.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kakavandi E, Shahbahrami R, Goudarzi H, Eslami G, Faghihloo E. Anoikis resistance and oncoviruses. J Cell Biochem (2018) 119(3):2484–91. doi: 10.1002/jcb.26363

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Li S, Chen Y, Zhang Y, Jiang X, Jiang Y, Qin X, et al. Shear stress promotes anoikis resistance of cancer cells via caveolin-1-dependent extrinsic and intrinsic apoptotic pathways. J Cell Physiol (2019) 234(4):3730–43. doi: 10.1002/jcp.27149

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhi Z, Ouyang Z, Ren Y, Cheng Y, Liu P, Wen Y, et al. Non-canonical phosphorylation of bmf by p38 MAPK promotes its apoptotic activity in anoikis. Cell Death Differ (2022) 29(2):323–36. doi: 10.1038/s41418-021-00855-3

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Jin L, Chun J, Pan C, Alesi GN, Li D, Magliocca KR, et al. Phosphorylation-mediated activation of LDHA promotes cancer cell invasion and tumour metastasis. Oncogene (2017) 36(27):3797–806. doi: 10.1038/onc.2017.6

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Jiang K, Yao G, Hu L, Yan Y, Liu J, Shi J, et al. MOB2 suppresses GBM cell migration and invasion via regulation of FAK/Akt and cAMP/PKA signaling. Cell Death Dis (2020) 11(4):230. doi: 10.1038/s41419-020-2381-8

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Kim H, Choi P, Kim T, Kim Y, Song BG, Park YT, et al. Ginsenosides Rk1 and Rg5 inhibit transforming growth factor-β1-induced epithelial-mesenchymal transition and suppress migration, invasion, anoikis resistance, and development of stem-like features in lung cancer. J Ginseng Res (2021) 45(1):134–48. doi: 10.1016/j.jgr.2020.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Du L, Han XG, Tu B, Wang MQ, Qiao H, Zhang SH, et al. CXCR1/Akt signaling activation induced by mesenchymal stem cell-derived IL-8 promotes osteosarcoma cell anoikis resistance and pulmonary metastasis. Cell Death Dis (2018) 9(7):714. doi: 10.1038/s41419-018-0745-0

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Yan H, Bu P. Non-coding RNA in cancer. Essays Biochem (2021) 65(4):625–39. doi: 10.1042/ebc20200032

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Park EG, Pyo SJ, Cui Y, Yoon SH, Nam JW. Tumor immune microenvironment lncRNAs. Brief Bioinform (2022) 23(1):bbab504. doi: 10.1093/bib/bbab504

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Smyth MJ, Ngiow SF, Ribas A, Teng MW. Combination cancer immunotherapies tailored to the tumour microenvironment. Nat Rev Clin Oncol (2016) 13(3):143–58. doi: 10.1038/nrclinonc.2015.209

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Kaymak I, Williams KS, Cantor JR, Jones RG. Immunometabolic interplay in the tumor microenvironment. Cancer Cell (2021) 39(1):28–37. doi: 10.1016/j.ccell.2020.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Dyson KA, Stover BD, Grippin A, Mendez-Gomez HR, Lagmay J, Mitchell DA, et al. Emerging trends in immunotherapy for pediatric sarcomas. J Hematol Oncol (2019) 12(1):78. doi: 10.1186/s13045-019-0756-z

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Chen C, Xie L, Ren T, Huang Y, Xu J, Guo W. Immunotherapy for osteosarcoma: Fundamental mechanism, rationale, and recent breakthroughs. Cancer Lett (2021) 500:1–10. doi: 10.1016/j.canlet.2020.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar