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.

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).

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:

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).

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.

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).

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.

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).

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.

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.

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 ).

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

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 ).

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.