Differential Gene and lncRNA Screening
22 metastasis and 66 primary OS samples were utilized in this study. To select differential genes and lncRNA, the Limma R package was adopted. As a result, 135 differential lncRNAs were identified, including 79 upregulated and 56 downregulated (Figure 1A). In addition, there were 171 differential genes (Figure 1B).
Construction of Prognosis Assessment Model
The results of the univariate cox regression analysis of 135 lncRNAs were used in the LASSO regression to identify robust markers. After the multivariate cox analysis, we identified seven prognosis-related lncRNAs (Figure 2A-C), which consisted of AL512422.1, AL008718.3, C5orf66-AS1, AL360182.2, CEBPA-DT, AC006033.2 and AL357507.1. All of these lncRNAs were taken part in prognosis model establishing, namely, y= 2.603*AL512422.1+3.735*AL008718.3+18.047* C5orf66-AS1+1.292*AL360182.2+0.475*CEBPA-DT+0.038* AC006033.2 +1.812*AL357507.1. The patients were divided into a high-risk and low-risk group based on the median risk score calculated by the above formula. The expression of seven prognosis-related lncRNAs in high and low-risk groups was depicted in Figure 2D. Moreover, the model performed stronger prediction power and all area under the curve (AUC) values were all >0.8 1 year: AUC = 0.92, 95% Cl = 0.83–1.01; 3 years: AUC = 0.87, 95% Cl = 0.79–0.96; 5 years: AUC = 0.86, 95% Cl = 0.76–0.96) (Figure 2D). According to KM survival curves (Figure 2D), the high-risk group had significantly poor overall survival time compared with the low-risk group (p<0.0001, HR=1.09, 95% CI = 1.06-1.13). In addition, the KM survival curves suggested that five lncRNAs, including AL512422.1, AL008718.3, C5orf66-AS1, AL360182.2 and AL357507.1, were related to poor overall survival. In contrast, CEBPA-DT and AC006033.2 were associated with better survival (Figure 3A-G).
Co‑expression network and enrichment analysis
To explore the mRNAs associated with predict model-related lncRNAs and their potential TF, we integrated the lncRNA and gene data to construct the co-expression network. Firstly, we selected the metastasis-related genes closely associated with seven prognosis-related lncRNAs above reported (R>0.3, p<0.001). As a result, a total of 6 prognosis-related lncRNAs and 24 metastasis-related genes were used to construct the co-expression network (Figure 4A). The AC006033.2 and AL512422.1 were the most important nodes. To further understand the functions of these genes, we analyzed the interaction between 24 metastasis-related genes and 13 TFs (Figure 4B, Table 1). Then, the GO analysis was carried out for functional enrichment of the nodes in this co-expression network (Figure 4C). The results showed several biological processes (BP) and cellular components (CC) pathways markedly enriched. From the BP perspective, the nodes were mainly enriched in “catecholamine metabolic process”, “lymphocyte proliferation”, “mononuclear cell proliferation”, “lymphocyte differentiation” and “leukocyte proliferation”. From the CC perspective, “integral component of synaptic membrane”, “intrinsic component of synaptic membrane” and “glutamatergic synapse” were highlighted.
GSVA Analysis of seven prognosis-related lncRNAs
We performed GSVA to determine the dynamics of biological processes and pathways for Hallmark gene sets based on seven prognosis-related lncRNAs. As Figure 5 showed, “G2M_CHECKPOINT”, “DNA_REPAIR” and “PI3K_AKT_MTOR_SIGNALING” were remarkably activated in the high AL512422.1 group, while the low AL512422.1 group was enriched in “WNT_BETA_CATENIN_SIGNALING” and “KRAS_SIGNALING_DN”. Upregulation of AL008718.3 activated “WNT_BETA_CATENIN_SIGNALING” and “TGF_BETA_SIGNALING”, but downregulation of AL008718.3 led to “HEDGEHOG_SIGNALING” and “PI3K_AKT_MTOR_SIGNALING”. In addition, AC006033.2 high level was significantly enriched for “HYPOXIA”, “APOPTOSIS” and “NOTCH_SIGNALING”, whereas low level triggered “DNA_REPAIR” and “G2M_CHECKPOINT” process. Notably, “GLYCOLYSIS” and “NOTCH_SIGNALING” were all enriched in AL357507.1 and AL360182.2 low groups. For the C5orf66-AS1 low group, “TGF_BETA_SIGNALING”, “HYPOXIA” and “NOTCH_SIGNALING” were enriched. The C5orf66-AS1 high group is related to “WNT_BETA_CATENIN_SIGNALING”. At last, we observed enrichment of “G2M_CHECKPOINT”, “DNA_REPAIR” and “PI3K_AKT_MTOR_SIGNALING” in the CEBPA-DT high group and “WNT_BETA_CATENIN_SIGNALING” in the low group. The top high-frequency enriched pathways were “NOTCH_SIGNALING”, “HYPOXIA”, “WNT_BETA_CATENIN_SIGNALING” and “HEDGEHOG_SIGNALING”.
Correlation of prognosis-related lncRNAs and infiltrating immunocyte fractions
To investigate the diversities between the primary and metastasis OS, 22 infiltrating immunocyte fractions were compared. As Figure 6A showed, the difference between the two groups was insignificant. Next, we analyzed the relationship between seven prognosis-related lncRNAs and 22 infiltrating immune cells (Figure 6B). The result suggested that AC006033.2 had the strongest relationship with resting dendritic cells, macrophages M0, macrophages M2, activated mast cells, activated neutrophils and memory T cells CD4. AL357507.1 had a close relationship with dendritic cells resting, macrophages M1, neutrophils and plasma cells. CEBPA-DT was associated with monocyte, T cells CD4 naïve, T cells CD8 and T cells regulatory (Tregs). C5orf66-AS1 and AL512422.1 were only related to dendritic cells activated and dendritic cells resting, respectively. On the contrary, for AL360182.2 and AL008718.3, there was no association with immune cells.
The mRNA expression of seven lncRNAs in OS cell lines
We evaluated the mRNA level of seven lncRNAs (AL512422.1, AL008718.3, C5orf66-AS1, AL360182.2, CEBPA-DT, AC006033.2 and AL357507.1) in five OS cell lines. MNNG cell line was a highly OS metastasis cell line, while the other four were weakly OS metastasis cell lines (Figure 7). The result revealed that AL512422.1 had a higher expression level in U2OS and SJSA-1 cell lines than MNNG (Figure 7A). A higher mRNA level of AL360182.2 and C5orf66-AS1 is observed in four weakly metastasis cell lines (Figure 7E-F). Instead, AL008718.3 mRNA level is significantly down-regulated in USO2, SAOS-2, SJSA-1 and HOS cell lines (Figure 7B). Additionally, the mRNA level of AL357507.1 was only up-regulated in the SAOS-2 cell line (Figure 7D). Interestingly, even though they were all weak metastasis cell lines, the expression pattern is completely different. AC006033.2 mRNA level is higher in U2OS, SAOS-2 and HOS cell lines, and lower in SJSA-1 cell lines (Figure 7C). A similar situation occurred with CEBPA-DT, the mRNA level is higher in HOS cell line and lower in SJSA-1 cell line (Figure 7G).