Identification of Immune-Related Prognostic Genes and LncRNAs Biomarkers Associated With Osteosarcoma Microenvironment

Osteosarcoma (OS) is the most common malignancy of the bone that occurs majorly in young people and adolescents. Although the survival of OS patients markedly improved by complete surgical resection and chemotherapy, the outcome is still poor in patients with recurrent and/or metastasized OS. Thus, identifying prognostic biomarkers that reflect the biological heterogeneity of OS could lead to better interventions for OS patients. Increasing studies have indicated the association between immune-related genes (IRGs) and cancer prognosis. In the present study, based on the data concerning OS obtained from TARGET (Therapeutically Applicable Research to Generate Effective Treatments) database, we constructed a classifier containing 12 immune-related (IR) long non-coding RNAs (lncRNAs) and 3 IRGs for predicting the prognosis of OS by using the least absolute shrinkage and selection operation Cox regression. Besides, based on the risk score calculated by the classifier, the samples were divided into high- and low-risk groups. We further investigated the tumor microenvironment of the OS samples by ESTIMATE and CIBERSORT algorithms between the two groups. Finally, we identified three small molecular drugs with potential therapeutic value for OS patients with high-risk score. Our results suggest that the IRGs and IR-lncRNAs–based classifier could be used as a reliable prognostic predictor for OS survival.

Introduction

Osteosarcoma (OS) is the most common malignancy of bone that occurs majorly in young people and adolescents (1), which accounts for ~4–5 per 1,000,000 per year (2). Surgical resection combined with chemotherapy has increased long-term survival rates to 60–70% for patients with the localized OS, but only 20–30% for patients with recurrent and/or metastasized OS (3, 4). Besides, the outcome of OS patients may be distinctly different even with the same stage. Thus, identifying prognostic biomarkers that reflect the biological heterogeneity of OS could lead to better interventions for patients.

Much attention has been paid to immune oncology for its impressive clinical benefits in a variety of malignancies. Immune-related genes (IRGs) and immune infiltrating cells have been considered as determining factors for regulating tumor development and progression (5, 6). The breakthrough of immunomodulatory therapies targeting the programmed death 1 (PD-1)/PD-1 ligand (PD-L1) signaling, a pathway crucial for impairing the immune system, has shown considerable success in multiple cancers by promoting antitumor immune function (7). Besides, studies have shown that both PD-1 and PD-L1 were significantly upregulated in OS patients and correlated with the prognosis (8, 9), and a recent study found that blockade of PD-1/PD-L1 signaling dramatically promoted the activity of cytotoxic T lymphocytes, inhibiting the tumor growth and increasing the survival rate in the mouse model of metastatic OS (10). Besides, tumor microenvironment (TME), where the tumor cells are located, is increasingly thought to play vital roles in tumor development and progression (11). Tumor microenvironment consists of extracellular matrix molecules, stromal cells, immune cells, and inflammatory mediators (12). As one of the major non-tumor cellular populations in the TME, infiltrating immune cells have been shown highly associated with responses to treatments and clinical outcomes of cancers. Tumor with high immune infiltration was associated with a better prognosis (13–16). Besides, bone has a highly specialized immune environment, and multiple immune signaling pathways play important roles in bone homeostasis (3). These evidences suggest that the application of immune-based prognostic biomarkers in OS is a potential. Furthermore, based on this, we can explore the underlying mechanisms and identify potential therapeutic drugs, so as to bring new insights into the improvement of the prognosis of OS patients.

Recent evidence suggests that long non-coding RNAs (lncRNAs) play important roles in regulating the development and activation of multiple immune cells through controlling the dynamic transcriptional programs (17) and involve carcinogenesis and metastasis (18). Immune-related (IR) lncRNAs have shown to act as biomarkers for predicting the risk of cancer patients of gliomas (19) and glioblastoma multiforme (20).

Several IRG-based signatures have been constructed to predict the risk of patients with different cancer types, such as lung cancer (21), glioblastoma (22), gastric cancer (23), and renal papillary cell carcinoma (24). As for OS, the prognostic values of IRGs and IR-lncRNA have still not been explored. In the present study, in an effort to assess the potential utility of IRGs and IR-lncRNAs in the prognosis of OS, we constructed a classifier containing 12 IR-lncRNAs and 3 IRGs for overall survival by using the least absolute shrinkage and selection operation (LASSO) Cox regression. Based on the risk score calculated by the classifier, the samples were divided into low- and high-risk groups. We further investigated the TME of the OS samples by Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression (ESTIMATE) data and Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithms. Finally, we explore the potential therapeutic small molecular drugs for the OS patients at high risk. Our results demonstrate that the IRGs and IR-lncRNAs–based classifier could be used as a reliable prognostic predictor for OS survival.

Materials and Methods

Data Source and Preprocessing

The bioinformatics analysis was conducted following the procedure presented in Figure 1. An RNA-seq data set and the corresponding clinical parameters of OS tissues (n = 88) were downloaded from TARGET (Therapeutically Applicable Research to Generate Effective Treatments) (https://ocg.cancer.gov/programs/target). The clinical characteristics of the 88 included samples are summarized in Table 1. The OS patients with complete outcome data and RNA-seq data were included in the subsequent analysis.

FIGURE 1

Figure 1. Flowchart of the study involved in construction of IRGs and IR-lncRNAs–based prognostic classifier. IRGs, immune related genes; TARGET, Therapeutically Applicable Research to Generate Effective Treatments; ROC, receiver operating characteristic; K-M, Kaplan–Meier; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; WGCNA, weighted gene coexpression network analysis; TME, tumor microenvironment; cMap, connectivity map; CIBERSORT, Cell Type Identification by Estimating Relative Subsets of RNA Transcripts.

TABLE 1

Table 1. Clinical characteristics of the osteosarcoma patients in TARGET database.

A total of 1,811 IRGs were obtained from Immport Shared Data (https://www.immport.org/shared/home).

Weighted Gene Co-expression Network Construction and Interesting Module Detection

Weighted gene co-expression network construction and module identification of all IRGs in the OS data set were performed following the protocols of weighted gene co-expression network analysis (WGCNA) (25), described previously (26). Briefly, every pairwise gene–gene relationship was calculated by a gene coexpression similarity in the first step. Then, a “soft” power adjacency function was utilized to construct the adjacency matrix and topological overlap matrix (TOM). “Gene modules,” groups of genes that have high topological overlap, were identified using hierarchical clustering with a dissimilarity measure (1-TOM).

The correlations between modules and clinical features were identified by Pearson correlation tests to identify clinically meaningful modules. The modules that exhibited a high correlation with prognostic features were selected as interesting modules to be further studied.

Identification of Prognostic IRGs and IR-lncRNAs

We conducted a univariate Cox regression for all IRGs in interesting modules and identified the genes with P < 0.05 as prognostic IRGs. Afterward, we conducted Pearson correlation tests between prognostic genes and all lncRNAs of the OS patients; correlation ≥0.6 was identified as IR-lncRNAs. Following this, we conducted a univariate Cox regression for all IR-lncRNAs and identified the lncRNAs with P < 0.05 as prognostic IR-lncRNAs.

Establishment of Prognostic Classifier

We applied the LASSO Cox regression analysis for all prognostic IRGs and IR-lncRNAs to select the most useful prognostic biomarkers and construct the survival-predicting classifier. The prognosis risk scores were calculated based on a formula as follows:

$\begin{array}{l}\text{Risk score}=\Sigma \text{ genes }(\text{or lncRNAs}){}_{\text{Cox coefficient}}\\ \times \text{ genes }(\text{or lncRNAs}){}_{\text{expression levels}}\end{array}$

Based on the cutoff of the median risk score, OS patients were divided into low- and high-risk groups. The predictive ability of the model for the training and validation cohort, which randomly split at a 1:1 ratio, as well as the total cohort, was evaluated by the Kaplan–Meier log-rank test. In addition, the application value of the model was tested by using receiver operating characteristic (ROC) curve analysis, and by univariate and multivariate Cox regression analysis.

Estimation of Immune Score

ESTIMATE was conducted to investigate the TME of OS and explore its correlation with IRGs and IR-lncRNAs–based classifier. ESTIMATE was designed to calculate scores for reflecting the levels of infiltrating immune cells and stromal cells within the TME based on the specific gene expression signatures of immune and stromal cells (27). Based on the cutoff of the median immune score, OS patients were divided into two groups. Besides, Kaplan–Meier method was also applied to assess the relationship between the immune score and the overall survival of OS patients.

Estimation of Immune Cell Infiltration

In order to further explore the relationship between the classifier and immune cell infiltration, the CIBERSORT algorithm was used to estimate the fraction of 22 immune cell types in the OS samples from gene expression data. Samples with a CIBERSORT output of P < 0.05 were considered to be eligible for further analysis. The Wilcoxon rank-sum test was used to identify the immune cells, which had significant differences in the proportion between low- and high-risk groups. Besides, Kaplan–Meier method was also applied to assess the relationship between the infiltrating of immune cells and the overall survival of OS patients.

Identification of Differentially Expressed Genes and Pathway Enrichment Analysis

The “edgeR” package of R was used to detect the differentially expressed gene (DEGs) between high- and low-risk samples. We set |log₂ fold change (FC)| ≥1 and P < 0.05 as the cutoff criteria. The volcano plot was drawn through the “gplots” package of R. Pathway enrichment analysis of DEGs, including KEGG pathway, Reactome pathway, and PANTHER pathway, was conducted by KOBAS 3.0 database (http://kobas.cbi.pku.edu.cn/anno_iden.php).

Identification of Potential Small Molecule Drugs

We submitted the DEGs of |log₂FC| ≥2 into the CMap (Connectivity map) (https://portals.broadinstitute.org/cmap), a database of small-molecule drugs, gene expression profiles, and diseases, which is based on the differential gene expression of human cells treated with small-molecule drugs. An enrichment score representing similarity is finally calculated. The positive connectivity score illustrates that the drug is capable of increasing the risk of death of OS patients. On the contrary, the negative link score indicated that the drug is able to reduce the risk of death. The drugs with negative connectivity score indicated potential therapeutic value. Two-dimensional diagrams of these candidate molecular drugs were obtained in Pubchem database (https://pubchem.ncbi.nlm.nih.gov/).

Statistical Analysis

All statistical analyses were performed by SPSS version 22.0, IBM Corp., Armonk, NY, United States and R version 3.6.1 software. The correlation between risk score and clinicopathological characteristics was analyzed by the χ² test. The unpaired t-test was used to estimate the statistical significance for normally distributed variables of the two groups. Wilcoxon rank-sum test was used to estimate the statistical significance for non–normally distributed variables of the two groups. All statistical tests were two-tailed, with a value of P < 0.05 considered statistically significant.

Results

Coexpression Network Construction and Interesting Module Detection

WGCNA was performed on 1,222 IRGs in the 88 OS samples. After removing one outlier sample, the connectivity between the genes in the gene network formed a scale-free network distribution when the soft-threshold power β was set to 8 (Figure 2A). Then, 10 coexpressed modules were identified and represented by different colors. The “gray” module was reserved for genes identified as not coexpressed (Figure 2B). The correlations between modules and clinical features, such as gender, race, age, EFS (event-free survival), overall survival, metastasis, progression, and death time were calculated. The red module was highly correlated with EFS (r = 0.33, P = 0.002), overall survival (r = 0.32, P = 0.002), and death time (r = 0.6, P = 6 × 10⁻⁶), and brown module was highly correlated with EFS (r = 0.34, P = 0.001), overall survival (r = 0.33, P = 0.002), and death time (r = 0.45, P = 6 × 10⁻⁵) (Figure 2C). Thus, the red and brown modules were selected as interesting modules to be studied in subsequent analyses.

FIGURE 2

Figure 2. WGCNA network and module detection. (A) Selection of the soft-thresholding powers. Power 8 was chosen because the fit index curve flattened out upon reaching a high value (>0.85). (B) Cluster dendrogram and module assignment for modules from WGCNA. The colored horizontal bar below the dendrogram represents the modules. (C) Correlation matrix for eigengene values and clinical features. Each cell includes the corresponding correlations and the p-values.

Identification of Prognostic IRGs and IR-lncRNAs

Eighty-six samples with complete survival data were included in the survival analysis. Univariate Cox regression analyses for all IRGs in red (n = 62) and brown (n = 180) modules were conducted (Supplementary Table 1) and identified 68 genes with P < 0.05 as prognostic IRGs (Supplementary Figure 1A). Afterward, 1,591 IR-lncRNAs were identified of correlation ≥0.6 with prognostic IRGs (Supplementary Table 2). Besides, 129 prognostic IR-lncRNAs were identified with P < 0.05 of univariate Cox regression analysis (Supplementary Figure 1B).

Establishment of Prognostic Classifier

LASSO Cox regression analysis was conducted to select the most useful prognostic biomarkers for constructing the prognostic-predicting classifier base on the training cohort (Figures 3A,B). A total of 12 lncRNAs (SNHG12, AL391421.1, AC117402.1, IL10RB-AS1, AL390038.1, AC083900.1, LINC01980, RUSC1-AS1, AC025822.1, AL133410.1, AL360182.2, and AL590764.1) and 3 mRNAs (IL7, SOCS1, and TMPRSS6) were identified as the most useful prognostic biomarkers (Table 2).

FIGURE 3

Figure 3. Construction of IRGs and IR-lncRNAs–based prognostic classifier. The results of the LASSO Cox regression suggested that all 15 mRNAs and lncRNAs were essential for the classifier (A,B). The expression levels of all 15 biomarkers of the classifier in high- and low-risk group from the training (C), validation (D), and total (E) cohorts. RS, risk score.

TABLE 2

Table 2. The IRGs and IR-lncRNAs in the prognostic classifier associated with OS in the TARGET data set.

The risk scores were calculated using the formula mentioned previously; patients in every cohort were further divided into high- and low-risk groups with the median risk score as the cutoff value. And the expression levels of every biomarker in different groups were analyzed (Figures 3C–E).

Correlation Between Classifiers and Clinicopathologic Characteristics

As shown in Table 3, all the clinical characteristics (age, gender, race, metastasis, and progression) showed no significant differences between the high- and low-risk groups in the training and validation cohort. However, metastasis showed significant difference between the two groups in the total cohort. Patients with metastasis were inclined to have a higher risk score.

TABLE 3

Table 3. Correlations between risk score of the immune-related genes and lncRNAs-based classifier with clinicopathological characteristics in training cohort, validation cohort and total cohort.

Prognostic Value of Classifiers for Assessing Overall Survival

As shown in Figures 4A–C, with the increase of risk score, the survival time of patients is decreased, and almost all the dead patients were enrolled in the high-risk group.

FIGURE 4

Figure 4. The prognostic value of IRGs and IR-lncRNAs–based classifier. The distribution of patients' risk score, survival state, and expression of all 15 biomarkers of the classifier in high- and low-risk group from the training (A), validation (B), and total (C) cohorts. Kaplan–Meier survival analysis of overall survival between high- and low-risk patients from the training (D), validation (E), and total (F) cohorts.

To further assess the prognostic value of the classifier, Kaplan–Meier test was conducted. As shown in Figures 4D–F, patients in high-risk group had significantly unfavorable prognosis.

Besides, the results of univariate Cox regression analysis in training, validation, and total cohorts further validated the prognostic value of classifier (Figures 5A–C). Moreover, multivariate analysis in the total cohort suggested that the classifier was an independent risk factor of survival for OS patients (Figure 5D).

FIGURE 5

Figure 5. The prognostic value of the classifier. Univariate and multivariate Cox regression analyses of the classifier with overall survival in the training (A), validation (B), and total (C,D) cohorts. The time-dependent ROC for 1-, 3-, and 5-years overall survival predictions for the classifier in training, validation, and total cohorts (E).

In addition, in the time-dependent ROC curve analysis, the areas under the curve for overall survival in the first, third, and fifth year were 1.009, 0.957, and 0.933, respectively, in the training cohort (Figure 5E), 0.945, 0.963, and 0.927, respectively, in the validation cohort; and 0.875, 0.956, and 0.927, respectively, in the total cohort. Moreover, the prediction capability of the classifier was superior to metastasis, which may be a major risk factor for tumor prognosis as reported by previous studies (Supplementary Figure 2).

The results above indicate that the IRGs and IR-lncRNAs–based classifier provided a useful prognostic tool with clinical value for appropriately categorizing patients with OS.

Patients With Low Risk Scores Correlated With High Immune Scores

As shown in Figure 6A, patients with low risk scores were related to high immune scores. Moreover, patients with high immune scores were correlated with better prognosis (Figure 6B).

FIGURE 6

Figure 6. (A) High-risk score correlated with low immune score. (B) Kaplan–Meier survival analysis of overall survival between high and low immune score. RS, risk score.

The Landscape of Immune Infiltration in OS

As shown in Figure 7A, we created a bar plot to demonstrate the proportion of 22 immune cells in each sample, which revealed that the five immune cells with the highest proportion in OS were M0 macrophages (38.6%), M2 macrophages (27.8%), T-cell CD4 memory resting (17.2%), mast cells resting (3.0%), and natural killer (NK) cells resting (2.9%). Then, we plotted the heat map of 22 immune cells in Figure 7B.

FIGURE 7

Figure 7. The composition (A) and heat map (B) of immune cells estimated by CIBERSORT algorithm in OSs. (C) The comparison of the fractions of immune cells between high- and low-risk group. (D) Kaplan–Meier survival analysis of overall survival between high and low level of infiltrating T-cell CD4 naive.

Additionally, Wilcoxon rank-sum test was used and revealed that the fractions of T-cell CD4 naive (P = 0.043), T-cell follicular helper (P = 0.049), dendritic cells resting (P = 0.049), and NK cells activated (P = 0.033) varied significantly between high- and low-risk-score patients (Figure 7C).

Besides, Kaplan–Meier analysis revealed that patients with low proportion of T-cell CD4 naive are associated with better overall survival (P = 0.05) (Figure 7D).

Screening for DEGs

A total of 1,135 DEGs, including 316 upregulated genes and 819 downregulated genes, were identified in the high-risk group, compared with the low-risk group (Figure 8A). We further performed pathway enrichment analysis for these DEGs. As shown in Figure 8B, the upregulated genes mainly enriched in the pathways of class A/1 (rhodopsin-like receptors), peptide ligand-binding receptors, GPCR ligand binding, GPCR downstream signaling, activation of C3 and C5, neuroactive ligand–receptor interaction, diseases of metabolism, signal transduction, inflammation mediated by chemokine and cytokine signaling pathway, and cell–cell communication. However, the downregulated genes mainly enriched in the pathways of transmembrane transport of small molecules, GPCR ligand binding, GPCR downstream signaling, class A/1 (rhodopsin-like receptors), neuroactive ligand–receptor interaction, starch and sucrose metabolism, retinol metabolism, drug metabolism–cytochrome P450, biological oxidations, and amino acid conjugation (Supplementary Figure 3).

FIGURE 8

Figure 8. Two-dimensional diagram of the three most significant drugs. (A) Thiamine, (B) harmalol, and (C) SC-19220.

Potential Small Molecule Drugs

We uploaded 404 DEGs of |log₂FC| ≥2, consisting of 300 downregulated genes and 104 upregulated genes, into the CMap network tool. Among these highly significant correlated molecules, thiamine, harmalol, and SC-19220 were most negatively correlated with high-risk OS patients (Figure 8). They all might have the potential therapeutic effects on OS patients with high risk.

Discussion

Osteosarcoma is the most common malignancy of bone and is characterized by highly aggressive and metastasis, which results in the very poor prognosis of patients (2). Thus, the identification of effective biomarkers for OS-specific prognoses is urgently needed to improve the management for patients. Taking into account the importance of the immune system in the progression of cancers and the highly specialized immune environment of bone (3, 28), it is essential to find out immune-related biomarker for the prognosis of OS patients, which may also play a significant role in immunotherapy.

In the present study, we constructed a prognostic classifier for OS by combining IRGs and IR-lncRNAs for the first time. A 12 IR-lncRNAs– and 3 IRGs-based classifiers for overall survival were constructed and validated to optimize the predictive ability of prognosis for OS patients. The results indicated that the classifier could successfully divide OS patients into high- and low-risk groups with significant differences in overall survival in the training cohorts. The prognostic value of the classifier was also confirmed in the validation cohort and the total cohort, indicating the repeatability and practicability of the IRGs- and IR-lncRNAs–based classifiers for the prognostic prediction for overall survival. Besides, the prediction capability of the classifier was superior to metastasis, which may be a major risk factor for tumor prognosis as reported by previous studies (29).

Among these 15 IR biomarkers, SNHG12, AL391421.1, AL390038.1, AC083900.1, RUSC1-AS1, AC025822.1, AL133410.1, and AL360182.2 were risk-associated, whereas AC117402.1, IL10RB-AS1, LINC01980, AL590764.1, IL7, SOCS1, and TMPRSS6 were protective (Table 1). Although some of the IR-lncRNAs in our classifier have not been functionally annotated and completely clarified, other biomarkers used in our classifiers have been explored. Previous studies showed that SNHG12 promoted tumorigenesis and metastasis in OS through upregulating NOCTH2 by sponging miR-195-5p (30). IL7 treatment promotes immune reconstitution significantly and improves the overall survival of pediatric sarcoma patients (31). SOCS1 acts as a cancer suppressor by promoting apoptosis and suppressing the metastasis of OS (32). Low expression of TMPRSS6 is related to the triple-negative and high grade of breast cancer (33). The upregulation of LINC01980 promotes tumor growth of esophageal squamous cell carcinoma (34). RUSC1-AS1 promotes the proliferation of breast cancer by inhibiting KLF2 and CDKN1A, which may serve as a potential hallmark for breast cancer (35). Given their strong relevance to prognosis, these genes should be explored in the future, especially in relation to OS.

Recently, many studies have demonstrated that tumor-infiltrating immune cells were associated with prognosis (36, 37). ESTIMATE algorithm is a simple method to predict the infiltration of immune cells by analyzing specific gene expression signature of immune cells and outputting immune scores (27). Previous ESTIMATE analyses have shown that immune cell infiltration is associated with prognosis in patients with various types of tumors (38, 39). In the present study, we found that the risk score based on the classifier negatively correlated with the immune score. Besides, patients with high immune scores have a favorable prognosis, indicating that immune cell infiltration of the TME could have a beneficial impact on prognosis. To further investigate the infiltration of immune cells, we conducted CIBERSORT analysis to illustrate the types of infiltrating cells. T-cell CD4 naive, T-cell follicular helper, dendritic cells resting, and NK cells activated were found significantly lower in the low-risk group. Besides, a low proportion of T-cell CD4 naive related to a better prognosis. Previous study revealed that tumor-infiltrating naive CD4 T cells are the important source of tumor-infiltrating regulatory T cells, which suppress the antitumor function of effector T cells and NK cells (40). Inhibiting the recruitment of T-cell CD4 naive into tumors reverses immunosuppression in breast cancer (41). Moreover, previous studies have shown that tumor-infiltrating T-cell follicular helper produced IL4 to suppress antitumor immunity by inducing myeloid cells to differentiate into M2 macrophages (42). Thus, infiltrating of T-cell CD4 naive and T-cell follicular helper may play important roles in the progression of OS, which will be well worth investigating.

Despite numerous attempts were done to improve the prognosis of OS, the outcome has remained unchanged over the past 3 decades (43). Herein, we identified three small molecules, thiamine, harmalol, and SC-19220, with potential therapeutic efficacy against OS. Thiamine (vitamin B₁) is a coenzyme for transketolase, pyruvate dehydrogenase, and α-ketoglutarate dehydrogenase complexes, which plays fundamental roles in various intracellular metabolism processes (44), as well as the regulation of immune system (45). The role of thiamine in immune responses has been demonstrated in the brain that thiamine plays significant anti-inflammation roles in inhibiting the expression of proinflammatory factors (cyclooxygenase-2, IL1, IL6, and TNF) and suppressing the CD40L-mediated immune and inflammatory responses (46). Current views on the role of thiamine in tumorigenesis are controversial (47). Some studies showed that thiamine was much higher in tumor tissues than in adjacent normal tissues (48), and a low dose of thiamine supplementations promoted cancer growth (49), suggesting that antithiamine is a potential way for cancer therapy. On the other hand, some studies showed that a high dose of thiamine reduced cell viability in breast cancer cells, but not in normal breast epithelial cells (50). Thus, the role of thiamine in OS is well worth investigating. Harmalol, a β-carboline alkaloid, presents in several plants such as Peganum harmala (51). Previous studies showed that harmalol treatment induced apoptosis of lung and liver cancer cells by activating caspase-8, caspase-3, and p53 (52, 53), indicating a potential antitumorous role of harmalol for OS. However, the role of harmalol on the immune system remains unclear. Prostaglandin E₂ (PGE₂) is a bioactive lipid that displays a wide array of biological effects associated with inflammation and cancer (54). Accumulation of PGE₂ in a cancer cell environment is a marker for the progression of many cancers (55). Blocking PGE₂ abrogates bladder cancer chemoresistance (56). SC-19220 is a prostaglandin E₂ antagonist, which showed potent anti-inflammation by suppression cytosolic phospholipase A₂ (57) and antitumor capacities by promoting tumor cell apoptosis through E-prostanoid 1 suppression (58). Collectively, thiamine, harmalol, and SC-19220 possess high clinical potential worthy of further investigation for the treatment of OS, especially through the mechanisms of modulating the immune system.

Inevitably, the present study has some innate limitations that need to be addressed. First, it was a retrospective study based on the publicly online database. Second, the cohort of the current study consisted of only 88 samples, and there is no cohort for validation from other databases. Thus, large-scale, multicenter studies are needed to confirm our results before the IRGs- and IR-lncRNAs–based classifier can be applied in the clinic.

Conclusion

In our study, we first identified and validated a classifier containing 12 IR-lncRNAs and 3 IRGs with independent prognostic significance for patients with OS. Moreover, our classifier can also provide novel clinical applications for immunotherapies and immune targets for OS. Besides, based on the classifiers, we identified three small molecular drugs with potential therapeutic value for OS treatment.

Data Availability Statement

Publicly available datasets were analyzed in this study, available from the TARGET database (https://ocg.cancer.gov/programs/target).

Author Contributions

TZ and YN: design, analysis and interpretation of data, drafting of the manuscript, and critical revision of the manuscript. HX and KC: statistical analysis. YZ: acquisition of data. XC, HL, and JW: critical revision of the manuscript for important intellectual content, administrative support, obtaining funding, and supervision. All authors: read and approved the final manuscript.

Funding

This work was supported by the Natural Science Foundation of China (Nos. 81570568, 81602419, and 81571075).

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.

The reviewer FL declared a shared affiliation, with no collaboration, with the authors to the handling editor at the time of the review.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.01109/full#supplementary-material

Supplementary Figure 1. Univariate Cox analysis for the IRGs (A) and IR-lncRNAs (B).

Supplementary Figure 2. The time-dependent ROC for 1-, 3-, and 5-years overall survival predictions for the classifier in comparison with clinical features in the training (A), validation (B), and total (C) cohorts.

Supplementary Figure 3. DEGs screen and pathway enrichment analysis. (A) Volcano plot of DEGs between high and low-risk groups in the total cohort. Pathway enrichment analysis of up- (B) and down- (C) regulated genes.

Supplementary Table 1. Gene list of red and brown modules.

Supplementary Table 2. List of Immune-related lncRNAs.

References

1. Hameed M, Mandelker D. Tumor syndromes predisposing to osteosarcoma. Adv Anat Pathol. (2018) 25:217–22. doi: 10.1097/PAP.0000000000000190

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ottaviani G, Jaffe N. The epidemiology of osteosarcoma. Cancer Treat Res. (2009) 152:3–13. doi: 10.1007/978-1-4419-0284-9_1

CrossRef Full Text | Google Scholar

3. Kansara M, Teng MW, Smyth MJ, Thomas DM. Translational biology of osteosarcoma. Nat Rev Cancer. (2014) 14:722–35. doi: 10.1038/nrc3838

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Bajpai J, Chandrasekharan A, Simha V, Talreja V, Karpe A, Pandey N, et al. Outcomes in treatment-naive patients with metastatic extremity osteosarcoma treated with OGS-12, a novel non-high-dose methotrexate-based, dose-dense combination chemotherapy, in a tertiary care cancer center. J Glob Oncol. (2018) 4:1–10. doi: 10.1200/JGO.17.00137

CrossRef Full Text | Google Scholar

5. Elinav E, Nowarski R, Thaiss CA, Hu B, Jin C, Flavell RA. Inflammation-induced cancer: crosstalk between tumours, immune cells and microorganisms. Nat Rev Cancer. (2013) 13:759–71. doi: 10.1038/nrc3611

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Gonzalez H, Hagerling C, Werb Z. Roles of the immune system in cancer: from tumor initiation to metastatic progression. Genes Dev. (2018) 32:1267–84. doi: 10.1101/gad.314617.118

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. (2012) 366:2443–54. doi: 10.1056/NEJMoa1200690

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Shen JK, Cote GM, Choy E, Yang P, Harmon D, Schwab J, et al. Programmed cell death ligand 1 expression in osteosarcoma. Cancer Immunol Res. (2014) 2:690–8. doi: 10.1158/2326-6066.CIR-13-0224

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yu F, Miao J. Significant association between cytotoxic T lymphocyte antigen 4 +49G>A polymorphism and risk of malignant bone tumors. Tumor Biology. (2013) 34:3371–5. doi: 10.1007/s13277-013-0908-7

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Lussier DM, O'Neill L, Nieves LM, McAfee MS, Holechek SA, Collins AW, et al. Enhanced T-cell immunity to osteosarcoma through antibody blockade of PD-1/PD-L1 interactions. J Immunothe. (2015) 38:96–106. doi: 10.1097/CJI.0000000000000065

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol. (2016) 27:1482–92. doi: 10.1093/annonc/mdw168

PubMed Abstract | CrossRef Full Text | Google Scholar

12. De Nola R, Menga A, Castegna A, Loizzi V, Ranieri G, Cicinelli E, et al. The crowded crosstalk between cancer cells and stromal microenvironment in gynecological malignancies: biological pathways and therapeutic implication. Int J Mol Sci. (2019) 20:2401. doi: 10.3390/ijms20102401

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Tas F, Erturk K. Tumor Infiltrating Lymphocytes (TILs) may be only an independent predictor of nodal involvement but not for recurrence and survival in cutaneous melanoma patients. Cancer Invest. (2017) 35:501–5. doi: 10.1080/07357907.2017.1351984

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Barnes TA, Amir E. HYPE or HOPE: the prognostic value of infiltrating immune cells in cancer. Br J Cancer. (2017) 117:451–60. doi: 10.1038/bjc.2017.220

CrossRef Full Text | Google Scholar

15. Karn T, Jiang T, Hatzis C, Sanger N, El-Balat A, Rody A, et al. Association between genomic metrics and immune infiltration in triple-negative breast cancer. JAMA Oncol. (2017) 3:1707–11. doi: 10.1001/jamaoncol.2017.2140

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Fridman WH, Zitvogel L, Sautes-Fridman C, Kroemer G. The immune contexture in cancer prognosis and treatment. Nat Rev Clin Oncol. (2017) 14:717–34. doi: 10.1038/nrclinonc.2017.101

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Atianand MK, Caffrey DR, Fitzgerald KA. Immunobiology of long non-coding RNAs. Annu Rev Immunol. (2017) 35:177–98. doi: 10.1146/annurev-immunol-041015-055459

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yang G, Lu X, Yuan L. LncRNA: a link between RNA and cancer. Biochimica et biophysica acta. (2014) 1839:1097–109. doi: 10.1016/j.bbagrm.2014.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Wang W, Zhao Z, Yang F, Wang H, Wu F, Liang T, et al. An immune-related lncRNA signature for patients with anaplastic gliomas. J Neurooncol. (2018) 136:263–71. doi: 10.1007/s11060-017-2667-6

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhou M, Zhang Z, Zhao H, Bao S, Cheng L, Sun J. An immune-related six-lncRNA signature to improve prognosis prediction of glioblastoma multiforme. Mol Neurobiol. (2018) 55:3684–97. doi: 10.1007/s12035-017-0572-9

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Song Q, Shang J, Yang Z, Zhang L, Zhang C, Chen J, et al. Identification of an immune signature predicting prognosis risk of patients in lung adenocarcinoma. J Transl Med. (2019) 17:70. doi: 10.1186/s12967-019-1824-4

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Cheng W, Ren X, Zhang C, Cai J, Liu Y, Han S, et al. Bioinformatic profiling identifies an immune-related risk signature for glioblastoma. Neurology. (2016) 86:2226–34. doi: 10.1212/WNL.0000000000002770

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Jiang B, Sun Q, Tong Y, Wang Y, Ma H, Xia X, et al. An immune-related gene signature predicts prognosis of gastric cancer. Medicine. (2019) 98:e16273. doi: 10.1097/MD.0000000000016273

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Wang Z, Song Q, Yang Z, Chen J, Shang J, Ju W. Construction of immune-related risk signature for renal papillary cell carcinoma. Cancer Med. (2019) 8:289–304. doi: 10.1002/cam4.1905

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhang T, Guo J, Gu J, Chen K, Wang Z, Li H, et al. KIAA0101 is a novel transcriptional target of FoxM1 and is involved in the regulation of hepatocellular carcinoma microvascular invasion by regulating epithelial-mesenchymal transition. J Cancer. (2019) 10:3501–16. doi: 10.7150/jca.29490

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med. (2015) 21:938–45. doi: 10.1038/nm.3909

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kempf-Bielack B, Bielack SS, Jurgens H, Branscheid D, Berdel WE, Exner GU, et al. Osteosarcoma relapse after combined modality therapy: an analysis of unselected patients in the Cooperative Osteosarcoma Study Group (COSS). J Clin Oncol. (2005) 23:559–68. doi: 10.1200/JCO.2005.04.063

CrossRef Full Text | Google Scholar

30. Zhou S, Yu L, Xiong M, Dai G. LncRNA SNHG12 promotes tumorigenesis and metastasis in osteosarcoma by upregulating Notch2 by sponging miR-195-5p. Biochem Biophys Res Commun. (2018) 495:1822–32. doi: 10.1016/j.bbrc.2017.12.047

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Merchant MS, Bernstein D, Amoako M, Baird K, Fleisher TA, Morre M, et al. Adjuvant immunotherapy to improve outcome in high-risk pediatric sarcomas. Clin Cancer Res. (2016) 22:3182–91. doi: 10.1158/1078-0432.CCR-15-2550

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Gong HL, Tao Y, Mao XZ, Song DY, You D, Ni JD. MicroRNA-29a suppresses the invasion and migration of osteosarcoma cells by regulating the SOCS1/NF-kappaB signalling pathway through negatively targeting DNMT3B. Int J Mol Med. (2019) 44:1219–32. doi: 10.3892/ijmm.2019.4287

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Tuhkanen H, Hartikainen JM, Soini Y, Velasco G, Sironen R, Nykopp TK, et al. Matriptase-2 gene (TMPRSS6) variants associate with breast cancer survival, and reduced expression is related to triple-negative breast cancer. Int J Cancer. (2013) 133:2334–40. doi: 10.1002/ijc.28254

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhang S, Liang Y, Wu Y, Chen X, Wang K, Li J, et al. Upregulation of a novel lncRNA LINC01980 promotes tumor growth of esophageal squamous cell carcinoma. Biochem Biophys Res Commun. (2019) 513:73–80. doi: 10.1016/j.bbrc.2019.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hu CC, Liang YW, Hu JL, Liu LF, Liang JW, Wang R. LncRNA RUSC1-AS1 promotes the proliferation of breast cancer cells by epigenetic silence of KLF2 and CDKN1A. Eur Rev Med Pharmacol Sci. (2019) 23:6602–11. doi: 10.26355/eurrev_201908_18548

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Jiang Y, Zhang Q, Hu Y, Li T, Yu J, Zhao L, et al. ImmunoScore signature: a prognostic and predictive tool in gastric cancer. Ann Surg. (2018) 267:504–13. doi: 10.1097/SLA.0000000000002116

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bense RD, Sotiriou C, Piccart-Gebhart MJ, Haanen J, van Vugt M, de Vries EGE, et al. Relevance of tumor-infiltrating immune cell composition and functionality for disease outcome in breast cancer. J Natl Cancer Inst. (2017) 109:djw192. doi: 10.1093/jnci/djw192

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Jia D, Li S, Li D, Xue H, Yang D, Liu Y. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging. (2018) 10:592–605. doi: 10.18632/aging.101415

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wang H, Wu X, Chen Y. Stromal-immune score-based gene signature: a prognosis stratification tool in gastric cancer. Front Oncol. (2019) 9:1212. doi: 10.3389/fonc.2019.01212

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zitvogel L, Tanchot C, Granier C, Tartour E. Following up tumor-specific regulatory T cells in cancer patients. Oncoimmunology. (2013) 2:e25444. doi: 10.4161/onci.25444

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Su S, Liao J, Liu J, Huang D, He C, Chen F, et al. Blocking the recruitment of naive CD4(+) T cells reverses immunosuppression in breast cancer. Cell Res. (2017) 27:461–82. doi: 10.1038/cr.2017.34

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Shirota H, Klinman DM, Ito SE, Ito H, Kubo M, Ishioka C. IL4 from T follicular helper cells downregulates antitumor immunity. Cancer Immunol Res. (2017) 5:61–71. doi: 10.1158/2326-6066.CIR-16-0113

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Bishop MW, Janeway KA, Gorlick R. Future directions in the treatment of osteosarcoma. Curr Opin Pediatr. (2016) 28:26–33. doi: 10.1097/MOP.0000000000000298

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Lu'o'ng KV, Nguyen LT. The role of thiamine in cancer: possible genetic and cellular signaling mechanisms. Cancer Genomics Proteomics. (2013) 10:169–85.

PubMed Abstract | Google Scholar

45. Spinas E, Saggini A, Kritas SK, Cerulli G, Caraffa A, Antinolfi P, et al. Crosstalk between vitamin B and immunity. J Biol Regul Homeost Agents. (2015) 29:283–8. doi: 10.6084/M9.FIGSHARE.3860199.V1

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Manzetti S, Zhang J, van der Spoel D. Thiamin function, metabolism, uptake, and transport. Biochemistry. (2014) 53:821–35. doi: 10.1021/bi401618y

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Boros LG, Brandes JL, Lee WN, Cascante M, Puigjaner J, Revesz E, et al. Thiamine supplementation to cancer patients: a double edged sword. Anticancer Res. (1998) 18(1b):595–602.

PubMed Abstract | Google Scholar

48. Baker H, Frank O, Chen T, Feingold S, DeAngelis B, Baker ER. Elevated vitamin levels in colon adenocarcinoma as compared with metastatic liver adenocarcinoma from colon primary and normal adjacent tissue. Cancer. (1981) 47:2883–6. doi: 10.1002/1097-0142(19810615)47:12 <2883::AID-CNCR2820471222>3.0.CO;2-I

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Comin-Anduix B, Boren J, Martinez S, Moro C, Centelles JJ, Trebukhina R, et al. The effect of thiamine supplementation on tumour proliferation. A metabolic control analysis study. Eur J Biochem. (2001) 268:4177–82. doi: 10.1046/j.1432-1327.2001.02329.x

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Liu X, Montissol S, Uber A, Ganley S, Grossestreuer AV, Berg K, et al. The effects of thiamine on breast cancer cells. Molecules. (2018) 23:1464. doi: 10.3390/molecules23061464

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Herraiz T, Gonzalez D, Ancin-Azpilicueta C, Aran VJ, Guillen H. beta-Carboline alkaloids in peganum harmala and inhibition of human monoamine oxidase (MAO). Food Chem Toxicol. (2010) 48:839–45. doi: 10.1016/j.fct.2009.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Abe A, Yamada H. Harmol induces apoptosis by caspase-8 activation independently of Fas/Fas ligand interaction in human lung carcinoma H596 cells. Anticancer Drugs. (2009) 20:373–81. doi: 10.1097/CAD.0b013e32832a2dd9

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Sarkar S, Bhattacharjee P, Bhadra K. DNA binding and apoptotic induction ability of harmalol in HepG2: biophysical and biochemical approaches. Chem Biol Interact. (2016) 258:142–52. doi: 10.1016/j.cbi.2016.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Nakanishi M, Rosenberg DW. Multifaceted roles of PGE2 in inflammation and cancer. Semin Immunopathol. (2013) 35:123–37. doi: 10.1007/s00281-012-0342-8

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Timoshenko AV, Xu G, Chakrabarti S, Lala PK, Chakraborty C. Role of prostaglandin E2 receptors in migration of murine and human breast cancer cells. Exp Cell Res. (2003) 289:265–74. doi: 10.1016/S0014-4827(03)00269-6

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Kurtova AV, Xiao J, Mo Q, Pazhanisamy S, Krasnow R, Lerner SP, et al. Blocking PGE2-induced tumour repopulation abrogates bladder cancer chemoresistance. Nature. (2015) 517:209–13. doi: 10.1038/nature14034

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lin CC, Lin WN, Wang WJ, Sun CC, Tung WH, Wang HH, et al. Functional coupling expression of COX-2 and cPLA2 induced by ATP in rat vascular smooth muscle cells: role of ERK1/2, p38 MAPK, and NF-kappaB. Cardiovasc Res. (2009) 82:522–31. doi: 10.1093/cvr/cvp069

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Chien CC, Ko CH, Shen SC, Yang LY, Chen YC. The role of COX-2/PGE2 in gossypol-induced apoptosis of colorectal carcinoma cells. J Cell Physiol. (2012) 227:3128–37. doi: 10.1002/jcp.23067

PubMed Abstract | CrossRef Full Text | Google Scholar