Bioinformatics analysis of lncRNAs in the occurrence and development of osteosarcoma

Abstract

Background: Osteosarcoma (OS) is a disease with high mortality in children and adolescents, and metastasis is one of its important clinical features. However, the molecular mechanism of OS occurrence is not completely clear. Thus, we screened potential biomarkers of OS and analyze their prognostic value.

Methods: The Cancer Genome Atlas (TCGA) datasets were used to analyze the differential lncRNAs in patients with OS of different immune score and the lncRNAs expressed by immune cells. Cox regression was used to develop the prognosis prediction model and specify the prognosis outcomes. Risk-proportional regression model was constructed, and the samples were divided into high and low groups based on the risk scores for the survival analysis. The areas under the receiver operating characteristic (ROC) curve were calculated and the risk-score model was verified. Finally, using 4 gene sets (comprising chemokines, immune checkpoint blockades, immune activity-related genes, and immune cells), and 4 analysis tools (CIBERSORT, TIMER, XCELL and MCP) to evaluated tumor immune infiltration.

Results: Twenty-nine long non-coding ribonucleic acids (lncRNAs) were obtained from the intersection of the screened lncRNAs. Caspase recruitment domain-containing protein 8-antisense RNA 1 (CARD8-AS1), lncRNA five prime to Xist (FTX), KAT8 regulatory NSL complex unit 1-antisense RNA 1 (KANSL1-AS1), Neuroplastin Intronic Transcript 1 (NPTN-IT1), oligodendrocyte maturation-associated long intervening non-coding RNA (OLMALINC) and RPARP Antisense RNA 1 (RPARP-AS1) were found to be correlated with survival. Univariate and multivariate regression analysis showed risk score [HR (hazard ratio) 3.5, P value 0.0043; HR 3.7, P value 0.0033] and metastasis (HR 4.7, P value 6.60E-05; HR 4.8, P value 8.36E-05) were the key factors of patients with OS. The areas under curves (AUCs) of the 1-, 3-, and 5-year ROC curves of the prognostic model were 0.715, 0.729, and 0.771. The low-risk patients tended to have a high abundance of immune cells.

Conclusions: This study showed that a risk score based on 6 lncRNAs has potential value in the prognosis of OS, and patients with low-risk scores have high immune cell infiltration and good prognosis. This study may enrich understandings of underlying mechanisms related to the occurrence and development of OS.

Keywords: Osteosarcoma (OS); biomarker; lncRNA; prognosis.

Figure 1

The volcano map of the DEGs between the samples with high and low immune scores were screened. Blue represents downregulated, red represents upregulated. The genes named in the figure are the lncRNAs expressed by the immune cells. The horizontal dotted line represents FDR0.05. The 2 vertical dashed lines represent –log2(1.5) and log2(1.5). DEGs, differentially expressed genes.

Figure 2

Prognostic value of the 6 lncRNAs in OS. (A) Survival curves for high- and low-risk patients. (B) ROC curves were drawn to analyze the predictive value of risk scores for the prognosis of OS. (C) The relationship between risk score and metastasis was analyzed. OS, osteosarcoma; ROC, receiver operating characteristics.

Figure 3

The DEGs in the two risk group samples were functionally analyzed. (A) The GO analysis of the biological processes. (B) The GO analysis of the cellular components. (C) The GO analysis of the molecular functions. Up gene (ratio): The ratio of the upregulated genes to total genes. If the value is 0.5, it means that the proportion of upregulated and downregulated genes is equal; if the value is >0.5, there are more upregulated genes. (D) KEGG analysis. The horizontal axis represents −log10 (P value); Color gradient fill according to −log10 (P value). The number on the right represents the number of differences; the colors are filled according to the number of downregulated genes. DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 4

A GSEA analysis was performed for the DEGs in the two risk group samples. (A) Immune-related gene sets were used for the analysis. (B) Gene sets from the KEGG database were used for the analysis. GSEA, Gene Set Enrichment Analysis; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; NES, Normalized Enrichment Score.

Figure 5

Chemokine and ICB gene sets were used to evaluate whether the model was associated with tumor immune invasion. (A) Boxplot of the chemokine gene concentration evaluation model. The horizontal axis is the name of the gene; the vertical axis is log2 expression. Each gene corresponds to 2 boxplots; that is, high risk and low risk, and are filled with blue and yellow, respectively. (B,D) Boxplot of the ICB gene concentration evaluation model. (C,E) The relationship between ICB gene expression level and risk score as presented in a heat map. Type stands for the risk score. Ns, not significant, *, P<0.05, **, P<0.01, ***, P<0.001, ****, P<0.0001. ICB, immune checkpoint blockade; TPM, transcript per million.

Figure 6

Immune activity-related and immune cell gene sets were used to evaluate whether the model was related to tumor immune invasion. (A) A boxplot of the evaluation model in the immune-activity-related gene set. (B) A boxplot of the evaluation model in the immune cell gene set. Ns, not significant, *, P<0.05, **, P<0.01, ***, P<0.001. TPM, transcript per million.

Figure 7

Immune cell infiltration analysis. CIBERSORT (A) and TIMER (B) were used to analyze the relationship between the risk score and immune infiltration results.

Figure 8

Immune cell infiltration analysis. XCELL (A) and MCP (B) were used to analyze the relationship between the risk score and immune infiltration results. (C) The correlations between the results of immune infiltration calculated by MCP and risk score were analyzed.