Epithelial to Mesenchymal Transition Relevant Subtypes with Distinct Prognosis and Responses to Chemo- or Immunotherapies in Osteosarcoma

Zhou, Yang; Li, Gai; Li, Hu; Lai, Fuchong; Duan, Pingguo; Cheng, Ming

doi:https://doi.org/10.1155/2022/1377565

Journal of Immunology Research

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusion Abbreviations Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Exploration of Novel Potential Biomarkers for Cancer Immunotherapy

View this Special Issue

Research Article | Open Access

Volume 2022 | Article ID 1377565 | https://doi.org/10.1155/2022/1377565

Epithelial to Mesenchymal Transition Relevant Subtypes with Distinct Prognosis and Responses to Chemo- or Immunotherapies in Osteosarcoma

Yang Zhou,¹Gai Li,¹Hu Li,¹Fuchong Lai,¹Pingguo Duan,¹and Ming Cheng¹

Academic Editor: Fu Wang

Received04 Apr 2022

Revised09 Jun 2022

Accepted13 Jun 2022

Published04 Jul 2022

Abstract

Objective. Currently, clinical classification of osteosarcoma cannot accurately predict the survival outcomes and responses to chemo- or immunotherapies. Our goal was to classify osteosarcoma patients into clinical/biological subtypes based on EMT molecules. Methods. This study retrospectively curated the RNA expression profiling of osteosarcoma patients from the TARGET and GSE21257 cohorts. Consensus clustering analyses were conducted in accordance with the expression profiling of prognostic EMT genes derived from univariate analyses. Immunological features were evaluated through immune cell infiltration, immune checkpoint expression, and activity of cancer immunity cycle. Drug sensitivity was estimated with the GDSC database. WGCNA approach was adopted to determine the EMT-derived genes. Following univariate analyses, a multivariate cox regression model was developed and externally verified. Predictive independency was evaluated with uni- and multivariate analyses. GSEA was presented to uncover relevant molecular mechanisms. Results. Prognostic EMT genes across osteosarcoma patients were stratified into distinct subtypes, namely, subtypes A and B. Patients in subtype B presented desirable prognosis, high immune activation, and enhanced sensitivity to cisplatin. Meanwhile, patients in subtype A were more sensitive to gemcitabine. In total, 86 EMT-derived hub genes were determined, and an EMT score was conducted for osteosarcoma prognosis. Following external verification, this EMT score was reliably and independently predictive of patients’ survival outcomes. Additionally, it was positively linked to steroid biosynthesis. Conclusion. Overall, our findings proposed EMT-relevant molecular subtypes and signatures for predicting prognosis and therapeutic responses, contributing to personalized treatment and clinical implication for osteosarcoma.

1. Introduction

Osteosarcoma represents the most prevalent primary bone sarcomas, which originates from mesenchymal stem cell populations, with an annual incidence of 3-5 cases per million individuals [1–3]. It most often affects children, adolescents, and young adults, with aberrant bone growth areas [4]. The difficulty in biology research is in relation to the complexity of the osteosarcoma genome as well as remarkable biological discrepancy among osteosarcoma subtypes [5]. Hence, an in-depth understanding of osteosarcoma biology highlights the heterogeneity as well as uncovers the molecular abnormality that defines patient subpopulations. The standard treatment of osteosarcoma comprises surgical resection and chemotherapy [6]. Patients are usually resistant to conventional chemotherapy; meanwhile, high-dose chemotherapeutic agents cause severe side effects [7]. Hence, to improve the survival duration of osteosarcoma patients has been proven to be challenging.

Epithelial-mesenchymal transition (EMT) is a key process by which epithelial cells acquire mesenchymal characteristics [8]. Growing evidences demonstrate that EMT exerts a crucial role in stemness [9], metabolic reprogramming [10], immune evasion [11], and therapeutic resistance [12] for cancer cells. Osteosarcoma cells that have experienced EMT process will lose the cellular polarity and acquire aggressive and metastatic features. This process has been regarded as a crucial event for osteosarcoma metastases [13]. Previously, Yiqi et al. proposed an EMT-relevant model of osteosarcoma as a prognostic indicator through integrated multicohorts [14], indicative of the crucial prognostic implication of EMT during osteosarcoma progression. Limited evidences demonstrate that chemotherapy (cisplatin, etc.) resistance contributes to EMT activation in osteosarcoma cells [15, 16]. Recent research has presented that EMT genes are remarkably linked to immunity in osteosarcoma [17]. Tumor-associated macrophages trigger EMT in osteosarcoma through activation of COX-2/STAT3 signaling [18]. Thus, extensive research of EMT on clinical outcomes and therapeutic responses of osteosarcoma is urgently required. In our research, integrated genomic analyses were conducted for characterizing distinct EMT-relevant subtypes with prognostic and prediction implications in osteosarcoma.

2. Materials and Methods

2.1. Acquisition of Gene Expression Profiling

This study acquired the TARGET data matrix (https://ocg.cancer.gov/programs/target/data-matrix) comprising gene expression data and clinical data of 84 osteosarcoma patients, as the discovery cohort. Genome-wide gene expression profiling and prognostic information of 53 osteosarcoma patients were curated from the GSE21257 cohort in the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/) [19], as the testing cohort. This cohort was on the basis of the platform of GPL10295 Illumina human-6 v2.0 expression beadchip. The gene expression series matrix files were directly retrieved, and the probe IDs were mapped to the gene symbols in accordance with the matched annotation files. Thereafter, expression values of multiple probes mapping to the same gene were averaged while probes mapping to multiple genes were eliminated. The gene set of EMT was curated from the Molecular Signatures Database (MSigDB), which was listed in Supplementary Table 1 [20].

2.2. Identification of Molecular Subtypes

Univariate analyses were presented for assessment of the association of EMT genes with osteosarcoma prognosis. Genes with were determined for the subsequent analyses. -means-based consensus clustering was carried out utilizing the ConsensusClusterPlus package on the basis of the transcriptome files of prognostic EMT genes [21]. Cumulative distribution function (CDF), delta area, and consensus matrix diagrams were conducted in accordance with default parameters. The number of clusters was determined in accordance with the following criteria: (i) the consistency within the clusters was relatively high; (ii) the coefficients of variation were relatively low; (iii) the area under the CDF curve was not remarkably elevated. The area under the CDF curves was utilized for defining the clustering number. Principal component analyses (PCA) were conducted for verifying the accuracy of this classification.

2.3. Gene Set Variation Analyses (GSVA)

Fifty hallmark pathways were acquired from the MSigDB, and the activity of above pathways was estimated with GSVA package in osteosarcoma specimens [22]. Limma package was adopted to compare the discrepancy in activity of hallmark pathways between subtypes [23].

2.4. Analyses of Immunological Features

The infiltration of tumor-infiltrating immune cells was estimated across osteosarcoma tissues utilizing the single-sample gene set enrichment analyses (ssGSEA) [24] and TIMER2 database [25]. The RNA expression of immune checkpoint molecules and HLA molecules was determined across osteosarcoma specimens. Cancer immunity cycle comprises release of cancer cell antigens (step 1), cancer antigen presentation (step 2), priming and activation (step 3), trafficking of immune cells to tumors (step 4), infiltration of immune cells into tumors (step 5), recognition of cancer cells by T cells (step 6), and killing of cancer cells (step 7) [26]. The activity of all steps was determined with ssGSEA approach in accordance with the transcriptome profiling [27]. The gene set of each step within the cancer immunity cycle was listed in Supplementary Table 2.

2.5. Estimation of Drug Sensitivity

Through adopting the pRRophetic algorithm [28], the half-maximal inhibitory concentration (IC50) values of cisplatin and gemcitabine were determined in osteosarcoma in accordance with the Genomics of Drug Sensitivity in Cancer (GDSC; http://www.cancerrxgene.org/) [29].

2.6. Weighted Gene Coexpression Network Analyses (WGNCA)

The data matrix of mRNA expression profiling in the TARGET cohort was analyzed utilizing the WGCNA package [30]. The first 25% genes with the largest variance were determined. For selecting a standard scale-free network, the sample hierarchical clustering approach was utilized for detecting and removing outlier specimens, followed by selection of the appropriate soft thresholding analyses. Thereafter, the adjacency matrix and topological overlap matrix (TOM) were established as well as the matched dissimilarity (1-TOM) was determined. Dynamic tree cutting methods were utilized for completing the gene tree and module clustering. Thereafter, the module characteristic genes were clustered as well as the highly similar modules were merged. The dissimilarity of the module eigengenes (ME) was determined with the moduleEigengenes function. The interactions of ME with EMT subtypes were analyzed with Pearson’s correlation. Module membership (MM) represents the relevance of the expression profile to ME, while gene significance (GS) represents the correlation between the expression profiles and clinical traits. In this study, MM was calculated with Pearson’s correlation of the expression profiling and ME. Meanwhile, GS was measured for evaluating the genes with EMT subtypes. Thereafter, the correlation between MM and gene signature was evaluated, and EMT-derived genes were determined.

2.7. Analysis of Hub Genes

Protein–protein interactions (PPIs) of genes in the lightcyan module were assessed through the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online tool (http://string-db.org/) [31]. Hub genes in this PPI network were determined utilizing the Molecular Complex Detection (MCODE), a plug-in of Cytoscape software [32].

2.8. Functional Enrichment Analyses

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were presented utilizing the clusterProfiler package [33]. The parameter utilized in this package was default as well as the threshold for identifying the GO functions and KEGG pathways was set as value < 0.05.

2.9. Gene Set Enrichment Analyses (GSEA)

For uncovering the difference in functional phenotypes between high- and low-risk subpopulations, GSEA was carried out. The “c5.all.v7.0.symbols.gmt” was utilized as a reference gene set. GSEA was implemented utilizing GSEA software. A nominal was regarded as a significant enrichment.

2.10. Statistical Analyses

Statistical analyses were implemented utilizing R software (version 4.0.2). Kaplan-Meier survival curves were drawn among osteosarcoma cases, and survival difference was determined with log-rank test. Receiver operator characteristic (ROC) curves were presented with timeROC package. Comparison between two subgroups was conducted with student’s or Wilcoxon test. Two-sided value < 0.05 was indicative of statistical significance.

3. Results

3.1. Characterization of Two EMT Molecular Subtypes in Osteosarcoma

Figure 1 illustrates the flowchart of this study. We curated 200 EMT genes from the MSigDB and analyzed their prognostic implication in osteosarcoma through univariate cox regression analyses. As a result, 40 EMT genes were remarkably linked to osteosarcoma prognosis (Figure 2(a); Table 1), in which 28 served as protective factors of osteosarcoma prognosis while the others served as risk factors. The consensus clustering of prognostic EMT genes was conducted to determine distinct EMT subtypes for osteosarcoma. For two categories, the area under the CDF curves began to stabilize (Figures 2(b) and 2(c)). Meanwhile, the consensus matrix was conducted for identifying the optimal number of subtypes. In Figure 2(d), the consensus matrix displayed a well-defined block structure when . Overall, two EMT molecular subtypes were conducted, namely, subtypes A and B. PCA results also fit into two clusters (Figure 2(e)). Survival analyses showed that subtype A presented poorer prognosis in comparison to subtype B (Figure 2(f)), indicative of the discrepancy in survival outcomes between subtypes.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2

Characterization of two EMT molecular subtypes across osteosarcoma patients in the TARGET cohort. (a) Forest plots depict the univariate cox regression analyses results of EMT genes that were significantly linked to osteosarcoma prognosis with . Red indicates risk factor while blue represents protective factor. (b) CDF curves when to 9. (c) The relative variations of the area under the CDF curves that from 2 to 9. (d) The consensus clustering matrix at . (e) PCA plots visualizing the discrepancy in two clusters in accordance with the expression profiling of prognostic EMT genes. (f) Kaplan-Meier survival analyses for patients in two clusters.

3.2. Immunological Features and Drug Sensitivity of Two EMT Molecular Subtypes

The activity of the 50 hallmark pathways was estimated across osteosarcoma tissues utilizing the GSVA algorithm. We noted that the remarkable activation of stromal pathways (EMT, Wnt β-catenin signaling, etc.), immune pathways (IL6 JAK STAT3 signaling, etc.), and metabolism pathways (bile acid metabolism, heme metabolism, xenobiotic metabolism, etc.) in EMT subtype B than subtype A (Figure 3(a) and Supplementary Table 3). In Figure 3(b), ssGSEA approach showed that EMT subtype B presented remarkably enhanced immune cell infiltrations containing central memory CD4 and CD8 T cells, effector memory CD4 T cells, memory B cells, regulatory T cells, type 1 and 2 helper cells, CD55bright natural killer cells, macrophages, MDSCs, natural killer cells, and natural killer T cells than subtype A. Meanwhile, TIMER2 approach showed the lower infiltrations of B cells and the higher infiltration of CD4 T cells compared with subtype A (Supplementary Figure 1). Additionally, we noted that immune checkpoint molecules ADORA2A, CD44, PDCD1LG2, TNFRSF14, and TNFRSF8 were prominently activated in subtype B than A (Figure 3(c)). Nevertheless, we did not investigate the significant discrepancy in HLA molecule expression between subtypes (Figure 3(d)). The activity of cancer immunity cycle was also estimated in each osteosarcoma specimen. Compared with subtype A, release of cancer cell antigen, B cell/monocyte/Th17 cell recruiting, recognition of tumor cells through T cells, and killer of tumor cells presented remarkably enhanced activity in subtype B (Figure 3(e)). Overall, EMT subtype B possessed the features of immune activation. Also, subtype B displayed increased sensitivity to cisplatin while subtype A was more sensitive to gemcitabine (Figure 3(f)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3

Immunological features and drug sensitivity of two EMT molecular subtypes. (a) Heat map depicts the activity of the 50 hallmark pathways in EMT subtypes A and B. (b) Boxplot shows the infiltrations of tumor-infiltrating immune populations in two subtypes. (c) Boxplot depicts the RNA expression of known immune checkpoint molecules in two subtypes. (d) Comparison of the RNA expression of HLA molecules in two subtypes. (e) Comparing the activities of all steps within cancer immunity cycle between subtypes. (f) Boxplot displays the sensitivity to cisplatin and gemcitabine in two subtypes. ; ; and .

3.3. Establishment of EMT Subtype-Relevant Coexpression Modules

Considering the sensitivity of WGCNA to batch effects, this study firstly preprocessed the data of osteosarcoma individuals in EMT molecular subtypes A and B from the TARGET cohort. The first 25% genes with the largest variance were determined for subsequent analyses. Thereafter, the hclust function was adopted for confirming the batch effect removal as well as investigating whether there was any outlier. In Figure 4(a), no outlier sample was found. Because of the premise of WGCNA approach required to hypothesize that the network met the scale-free criteria, this study further screened the optimal soft thresholding value for making a coexpression network was more subject to a scale-free network. Through calculation of the scale-free topology fitting index, value was up to 0.85 (Figure 4(b)), which validated the feasibility of WGCNA. A hierarchical clustering tree was retrieved through adopting hierarchical clustering analyses. Thereafter, in accordance with the dynamic tree cut approach, 29 distinct coexpression modules were assigned (Figure 4(c)). Figure 4(d) demonstrated that the lightcyan module presented the strongest correlation to EMT molecular subtypes. The genes in this module deserved more exploration.

(a)

(b)

(c)

(d)

Figure 4

Establishment of EMT subtype-relevant coexpression modules in the TARGET cohort. (a) Hierarchical clustering dendrogram of osteosarcoma individuals as well as the matched EMT molecular subtypes. (b) Analyses of network topology and average network connectivity under diverse soft-threshold powers (β). The red line is indicative of a correlation coefficient of 0.85. (c) Clustering dendrogram with dissimilarity in accordance with topological overlapping, along with designed module colors. (d) The interaction network between modules and EMT molecular subtypes. Each column is indicative of EMT molecular subtype while each row is indicative of a coexpression module. The number in the rectangle displays the correlation coefficient and value.

3.4. Exploration of EMT-Derived Genes and Their Relevant Biological Implication

Our further analyses uncovered that the genes in the lightcyan module were remarkably linked to two EMT molecular subtypes (Figures 5(a) and 5(b)), indicative of the genes in the lightcyan module as EMT-derived genes. Through MCODE analyses, we determined 86 EMT-derived hub genes, as depicted in Figure 5(c). Their relevant biological implication was explored in-depth. In Figure 5(d), the EMT-derived hub genes were remarkably linked to endoplasmic reticulum. Additionally, they were significantly enriched in ribosome, oxidative phosphorylation, chemical carcinogenesis-reactive oxygen species, thermogenesis, and HIF-1 signaling pathway (Figure 5(e)).

(a)

(b)

(c)

(d)

(e)

3.5. Development of an EMT-Derived Prognostic Model for Osteosarcoma

Our univariate cox regression analyses were indicative that 33 EMT-derived genes displayed significant interactions with osteosarcoma prognosis (Table 2). Above genes were input the multivariate cox regression model. In accordance with the expression and coefficient of candidate genes (Table 3), we developed an EMT-derived prognostic model for osteosarcoma, comprising RPS9, RPS23, EIF4A1, RPL12, RPL36, RPL37A, RPL34, EEF1B2, RPS8, RPS28, RPL10, RPS24, RPL35A, RPL11, RPL21, RPS27A, RPS12, and RPL13A. We noted the remarkable discrepancy in their expressions in high- than low-risk subpopulations (Figure 6(a)). Figure 6(b) displayed the EMT score distribution across osteosarcoma patients. As EMT score elevated, dead cases were gradually increased (Figure 6(c)). Survival analyses demonstrated that high-risk cases were indicative of more unfavorable clinical outcomes in comparison to low-risk cases (Figure 6(d)). ROC curves were presented to assess the predictive performance of EMT score in osteosarcoma prognosis. In Figure 6(e), area under the curves (AUCs) at one-, three-, and five-year survival were separately 0.723, 0.858, and 0.860, proving that EMT score was capable of prediction of osteosarcoma prognosis. Compared with previously published two EMT models from Yiqi et al. [14] and Peng et al. [17], our EMT-derived prognostic model had higher C-index (0.828), indicating the prediction superiority of this model (Figure 6(f)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6

Development of EMT score for osteosarcoma prognosis in the TARGET cohort. (a) Heat map depicts the expression of each gene from EMT score in each osteosarcoma specimen. Red is indicative of upregulation; and blue is indicative of downregulation. (b) The distribution of risk score among osteosarcoma individuals. Vertical dotted line is indicative of the median value of EMT score. (c) The distribution of survival status of high- and low-risk subgroups. (d) Kaplan-Meier survival analyses for two subgroups. (e) ROC curves at one-, three-, and five-year survival in accordance with EMT score. (f) C-index of our EMT model and previously published two EMT models from Zhang et al. and Peng et al.

3.6. External Verification of the Prognostic Value of EMT Score

The prognostic value of EMT score was further verified in the GSE21257 cohort. Figure 7(a) depicted the expression patterns of each gene from EMT score across osteosarcoma specimens. We also visualized the distribution of EMT score in the GSE21257 cohort (Figure 7(b)). In accordance with the median value of EMT score, we stratified osteosarcoma individuals into high- as well as low-risk subpopulations. Additionally, we noted that there were relatively increased dead cases in high-risk subgroup (Figure 7(c)). As expected, high EMT score was remarkably linked to more unfavorable survival outcomes (Figure 7(d)). The AUC at five-year survival was 0.818, proving the excellent predictive performance of EMT score (Figure 7(e)).

(a)

(b)

(c)

(d)

(e)

Figure 7

External verification of the prognostic value of EMT score in the GSE21257 cohort. (a) Heat map depicts the expression of each gene from EMT score in each osteosarcoma specimen. Red is indicative of upregulation while blue is indicative of downregulation. (b) The distribution of risk score among osteosarcoma individuals. Vertical dotted line is indicative of the mean value of EMT score. (c) The distribution of survival status of high- and low-risk subpopulations. (d) Kaplan-Meier survival analyses for two subpopulations. (e) ROC curves at five-year survival in accordance with EMT score.

3.7. Genes in EMT Score as Prognostic Indicators of Osteosarcoma

Further analyses were conducted for assessing the survival implication of each gene in EMT score in the TARGET cohort. As a result, high expression of EEF1B2, EIF4A1, RPL10, RPL11, RPL12, RPL13A, RPL21, RPL34, RPL35A, RPL36, RPL37A, RPS8, RPS9, RPS12, RPS23, RPS24, RPS27A, and RPS28 was prominently linked to more undesirable survival outcomes in comparison to their low expression (Figures 8(a)–8(r)).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

3.8. EMT Score Is Independent of Clinicopathological Indicators for Osteosarcoma Prognosis

In the TARGET cohort, our univariate cox regression analyses displayed that this EMT score was prominently linked to osteosarcoma survival outcomes (Figure 9(a)). Further multivariate analyses demonstrated that the EMT score was independently predictive of patients’ prognosis (Figure 9(b)). For uncovering the underlying biological phenotypes between subpopulations, this research carried out GSEA. As a result, steroid biosynthesis was remarkably activated in high-risk subpopulation (Figure 9(c)), and hypertrophic cardiomyopathy displayed prominent activation in low-risk subpopulation (Figure 9(d)).

(a)

(b)

(c)

(d)

4. Discussion

Osteosarcoma presents diverse clinical courses and biological heterogeneity [34]. Hence, to reliably predict prognosis and therapeutic responses, it is critical for comprehensively investigating the molecular mechanisms. Our study determined forty prognostic EMT genes in osteosarcoma. Through consensus clustering approach, two EMT molecular subtypes were conducted. Especially, EMT subtype A presented poorer prognosis in comparison to subtype B, indicating that there was a remarkable discrepancy in survival outcomes between subtypes. Recently, immunotherapy is a promising treatment option against osteosarcoma, and it is crucial to improve the comprehending about the immune responses [3]. We evaluated the immunological features from multiple perspectives. GSVA demonstrated the activation of immune pathways (like IL6 JAK STAT3 signaling) in subtype B. Our ssGSEA revealed that subtype B displayed remarkably enhanced immune cell infiltrations containing central memory CD4 and CD8 T cells, effector memory CD4 T cells, memory B cells, regulatory T cells, type 1 and 2 helper cells, CD55bright natural killer cells, macrophages, MDSCs, natural killer cells, and natural killer T cells. In addition, subtype B had the features of increased expression of immune checkpoint molecules ADORA2A, CD44, PDCD1LG2, TNFRSF14, and TNFRSF8. Among steps within cancer immunity cycle, release of cancer cell antigen, B cell/monocyte/Th17 cell recruiting, recognition of tumor cells through T cells, and killer of tumor cells presented remarkably enhanced activity in subtype B. Hence, EMT subtype B was characterized by immune activation. Drug resistance severely hinders the improvement of survival rate of osteosarcoma patients [35]. Cisplatin, an alkylating drug, can form irreversible covalent bonds with DNA, which causes DNA strands to cross-link and break and missense mutation [36]. It is widely applied for osteosarcoma chemotherapy as well as cisplatin resistance is frequent across osteosarcoma individuals [37]. Our data demonstrated that EMT subtype B presented enhanced sensitivity to cisplatin. Gemcitabine represents the second cytidine analog developed following cytosine arabinoside [38]. Intriguingly, subtype A was more sensitive to gemcitabine.

Through WGCNA approach, we constructed 29 EMT subtype-relevant coexpression modules in the TARGET cohort. Especially, lightcyan module presented the strongest interactions with EMT molecular subtypes. Further, MCODE analyses determined 86 EMT-derived hub genes. Functional enrichment analyses unraveled the remarkable interactions of these EMT-derived hub genes with endoplasmic reticulum (ER), ribosome, oxidative phosphorylation, chemical carcinogenesis-reactive oxygen species, thermogenesis, and HIF-1 signaling pathway. For instance, the ER represents the central intracellular organelle of diverse cellular functions and endoplasmic reticulum stress response participates in osteosarcoma pathogenesis [39]. SENP1/HIF-1alpha regulates hypoxia-mediated EMT in osteosarcoma cells [40].

This study conducted an EMT-derived prognostic model for osteosarcoma in the TARGET cohort, comprising EEF1B2, EIF4A1, RPL10, RPL11, RPL12, RPL13A, RPL21, RPL34, RPL35A, RPL36, RPL37A, RPS8, RPS9, RPS12, RPS23, RPS24, RPS27A, and RPS28. ROC curves proved that the EMT score was capable of accurately predicting osteosarcoma individuals’ clinical outcomes. Additionally, external verification proved the feasibility of this model in the GSE21257 cohort. Each gene in the EMT score was linked to undesirable prognosis of osteosarcoma individuals. Previously, mTOR inhibitor blunts the p53 response to nucleolar stress through modulating RPL11 and MDM2 expressions [41]. RPL34 upregulation is indicative of undesirable clinical outcomes in osteosarcoma as well as silencing RPL34 weakens osteosarcoma proliferation [42]. Suppressing RPS9 blunts osteosarcoma cell growth via inactivating MAPK signaling [43]. Nevertheless, the functions of most genes in osteosarcoma progression are lack of experimental evidences. Multivariate analyses suggested the independency of the EMT score in prediction of osteosarcoma outcomes. Further, GSEA results presented that steroid biosynthesis was remarkably activated in high-risk osteosarcoma individuals, indicating that steroid biosynthesis might affect osteosarcoma progression.

A few limitations of our study should be taken into consideration. First, further analyses are required for elucidating the molecular mechanisms underlying the impact of genes in the EMT-derived model on osteosarcoma through in-depth experiments. Second, more osteosarcoma patients can produce more convincing and accurate findings. Hence, abundant specimens will be included for improving our conclusions as well as more reliably illustrating the underlying mechanisms by which genes in the EMT-derived model affect osteosarcoma progression.

5. Conclusion

Collectively, this study adopted consensus clustering analyses for determining two EMT molecular subtypes in the TARGET cohort as well as uncovered the discrepancy in survival outcomes, immunological features, and drug sensitivity between subtypes. Through system biology-based WGCNA approach, we determined EMT-derived genes and conducted an EMT score for predicting osteosarcoma prognosis. To our knowledge, we firstly characterized the EMT-relevant subtypes for osteosarcoma. Additionally, the EMT score we proposed was externally verified in the external cohort. Overall, our findings might contribute to personalized treatment and be of much clinical implication for osteosarcoma.

Abbreviations

EMT:	Epithelial-mesenchymal transition
GEO:	Gene Expression Omnibus
MSigDB:	Molecular Signatures Database
CDF:	Cumulative distribution function
PCA:	Principal component analyses
GSVA:	Gene set variation analyses
ssGSEA:	Single-sample gene set enrichment analysis
GDSC:	Genomics of drug sensitivity in cancer
WGNCA:	Weighted gene coexpression network analyses
TOM:	Topological overlap matrix
ME:	Module eigengenes
MM:	Module membership
GS:	Gene significance
PPI:	Protein–protein interaction
STRING:	Search Tool for the Retrieval of Interacting Genes/Proteins
MCODE:	Molecular Complex Detection
GO:	Gene Ontology
KEGG:	Kyoto Encyclopedia of Genes and Genomes
GSEA:	Gene set enrichment analyses
ROC:	Receiver operator characteristic
AUCs:	Area under the curves.

Data Availability

The data used to support the findings of this study are included within the supplementary information files.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

Yang Zhou and Gai Li contributed equally to this work.

Acknowledgments

This work was funded by National Natural Science Foundation of China (81760393 and 81401790) and Natural Science Foundation of Jiangxi Province (Grant no. 20202BABL206034).

Supplementary Materials

Supplementary 1. Supplementary Figure 1: the differences in immune cell infiltrations between two EMT subtypes using TIMER2 approach.

Supplementary 2. Supplementary Table 1: the information of EMT genes from the Molecular Signatures Database.

Supplementary 3. Supplementary Table 2: the gene set of each step within the cancer immunity cycle.

Supplementary 4. Supplementary Table 3: comparison of the activity of hallmark pathways in two EMT molecular subtypes.

References

C. Du, M. Zhou, F. Jia et al., “D-arginine-loaded metal-organic frameworks nanoparticles sensitize osteosarcoma to radiotherapy,” Biomaterials, vol. 269, article 120642, 2021.
View at: Publisher Site | Google Scholar
Y. Yuan, G. Yan, M. He et al., “ALKBH5 suppresses tumor progression via an m⁶A-dependent epigenetic silencing of pre-miR-181b-1/YAP signaling axis in osteosarcoma,” Cell Death & Disease, vol. 12, no. 1, p. 60, 2021.
View at: Publisher Site | Google Scholar
C. Chen, L. Xie, T. Ren, Y. Huang, J. Xu, and W. Guo, “Immunotherapy for osteosarcoma: fundamental mechanism, rationale, and recent breakthroughs,” Cancer Letters, vol. 500, pp. 1–10, 2021.
View at: Publisher Site | Google Scholar
R. Rathore and B. A. van Tine, “Pathogenesis and current treatment of osteosarcoma: perspectives for future therapies,” Journal of Clinical Medicine, vol. 10, no. 6, article 1182, 2021.
View at: Publisher Site | Google Scholar
I. Fernandes, C. Melo-Alvim, R. Lopes-Brás, M. Esperança-Martins, and L. Costa, “Osteosarcoma pathogenesis leads the way to new target treatments,” International Journal of Molecular Sciences, vol. 22, no. 2, p. 813, 2021.
View at: Publisher Site | Google Scholar
J. Gill and R. Gorlick, “Advancing therapy for osteosarcoma,” Nature Reviews. Clinical Oncology, vol. 18, no. 10, pp. 609–624, 2021.
View at: Publisher Site | Google Scholar
C. Cheng, Q. Ding, Z. Zhang et al., “PTBP1 modulates osteosarcoma chemoresistance to cisplatin by regulating the expression of the copper transporter SLC31A1,” Journal of Cellular and Molecular Medicine, vol. 24, no. 9, pp. 5274–5289, 2020.
View at: Publisher Site | Google Scholar
X. Qiao, X. Niu, J. Shi et al., “Wnt5a regulates ameloblastoma cell migration by modulating mitochondrial and cytoskeletal dynamics,” Journal of Cancer, vol. 11, no. 18, pp. 5490–5502, 2020.
View at: Publisher Site | Google Scholar
Q. Tang, J. Chen, Z. Di et al., “TM4SF1 promotes EMT and cancer stemness via the Wnt/β-catenin/SOX2 pathway in colorectal cancer,” Journal of Experimental & Clinical Cancer Research, vol. 39, no. 1, p. 232, 2020.
View at: Publisher Site | Google Scholar
W. Hua, P. Ten Dijke, S. Kostidis, M. Giera, and M. Hornsveld, “TGFβ-induced metabolic reprogramming during epithelial-to-mesenchymal transition in cancer,” Cellular and Molecular Life Sciences, vol. 77, no. 11, pp. 2103–2123, 2020.
View at: Publisher Site | Google Scholar
R. Ren, Y. Du, X. Niu, and R. Zang, “ZFPM2-AS1 transcriptionally mediated by STAT1 regulates thyroid cancer cell growth, migration and invasion via miR-515-5p/TUSC3,” Journal of Cancer, vol. 12, no. 11, pp. 3393–3406, 2021.
View at: Publisher Site | Google Scholar
X. Zhu, L. Chen, L. Liu, and X. Niu, “EMT-mediated acquired EGFR-TKI resistance in NSCLC: mechanisms and strategies,” Frontiers in Oncology, vol. 9, p. 1044, 2019.
View at: Publisher Site | Google Scholar
X. Jiang, Z. Zhang, C. Song et al., “Glaucocalyxin A reverses EMT and TGF-β1-induced EMT by inhibiting TGF-β1/Smad2/3 signaling pathway in osteosarcoma,” Chemico-Biological Interactions, vol. 307, pp. 158–166, 2019.
View at: Publisher Site | Google Scholar
Z. Yiqi, L. Ziyun, F. Qin, W. Xingli, and Y. Liyu, “Identification of 9-gene epithelial-mesenchymal transition related signature of osteosarcoma by integrating multi cohorts,” Technology in Cancer Research & Treatment, vol. 19, 2020.
View at: Publisher Site | Google Scholar
E. Qiu, Y. Gao, B. Zhang, T. Xia, Z. Zhang, and G. Shang, “Upregulation of cell division cycle 20 in cisplatin resistance-induced epithelial-mesenchymal transition in osteosarcoma cells,” American Journal of Translational Research, vol. 12, no. 4, pp. 1309–1318, 2020.
View at: Google Scholar
S. Fang, L. Yu, H. Mei et al., “Cisplatin promotes mesenchymal-like characteristics in osteosarcoma through snail,” Oncology Letters, vol. 12, no. 6, pp. 5007–5014, 2016.
View at: Publisher Site | Google Scholar
Y. X. Peng, B. Yu, H. Qin, L. Xue, Y. J. Liang, and Z. X. Quan, “EMT-related gene expression is positively correlated with immunity and may be derived from stromal cells in osteosarcoma,” PeerJ, vol. 8, article e8489, 2020.
View at: Publisher Site | Google Scholar
Y. Han, W. Guo, T. Ren et al., “Tumor-associated macrophages promote lung metastasis and induce epithelial- mesenchymal transition in osteosarcoma by activating the COX-2/STAT3 axis,” Cancer Letters, vol. 440-441, pp. 116–125, 2019.
View at: Publisher Site | Google Scholar
E. P. Buddingh, M. L. Kuijjer, R. A. Duim et al., “Tumor-infiltrating macrophages are associated with metastasis suppression in high-grade osteosarcoma: a rationale for treatment with macrophage activating agents,” Clinical Cancer Research, vol. 17, no. 8, pp. 2110–2119, 2011.
View at: Publisher Site | Google Scholar
A. Liberzon, C. Birger, H. Thorvaldsdóttir, M. Ghandi, J. P. Mesirov, and P. Tamayo, “The molecular signatures database hallmark gene set collection,” Cell Systems, vol. 1, no. 6, pp. 417–425, 2015.
View at: Publisher Site | Google Scholar
M. D. Wilkerson and D. N. Hayes, “ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking,” Bioinformatics, vol. 26, no. 12, pp. 1572-1573, 2010.
View at: Publisher Site | Google Scholar
S. Hänzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and RNA-seq data,” BMC Bioinformatics, vol. 14, no. 1, p. 7, 2013.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015.
View at: Publisher Site | Google Scholar
G. Bindea, B. Mlecnik, M. Tosolini et al., “Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer,” Immunity, vol. 39, no. 4, pp. 782–795, 2013.
View at: Publisher Site | Google Scholar
T. Li, J. Fu, Z. Zeng et al., “TIMER2.0 for analysis of tumor-infiltrating immune cells,” Nucleic Acids Research, vol. 48, no. W1, pp. W509–w514, 2020.
View at: Publisher Site | Google Scholar
L. Chen, X. Niu, X. Qiao et al., “Characterization of interplay between autophagy and ferroptosis and their synergistical roles on manipulating immunological tumor microenvironment in squamous cell carcinomas,” Frontiers in Immunology, vol. 12, article 739039, 2022.
View at: Publisher Site | Google Scholar
L. Xu, C. Deng, B. Pang et al., “TIP: a web server for resolving tumor immunophenotype profiling,” Cancer Research, vol. 78, no. 23, pp. 6575–6580, 2018.
View at: Publisher Site | Google Scholar
P. Geeleher, N. Cox, and R. S. Huang, “pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels,” PLoS One, vol. 9, no. 9, article e107468, 2014.
View at: Publisher Site | Google Scholar
W. Yang, J. Soares, P. Greninger et al., “Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells,” Nucleic Acids Research, vol. 41, pp. D955–D961, 2013.
View at: Publisher Site | Google Scholar
P. Langfelder and S. Horvath, “WGCNA: an R package for weighted correlation network analysis,” BMC Bioinformatics, vol. 9, no. 1, p. 559, 2008.
View at: Publisher Site | Google Scholar
D. Szklarczyk, J. H. Morris, H. Cook et al., “The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible,” Nucleic Acids Research, vol. 45, no. D1, pp. D362–d368, 2017.
View at: Publisher Site | Google Scholar
N. T. Doncheva, J. H. Morris, J. Gorodkin, and L. J. Jensen, “Cytoscape StringApp: network analysis and visualization of proteomics data,” Journal of Proteome Research, vol. 18, no. 2, pp. 623–632, 2019.
View at: Publisher Site | Google Scholar
G. Yu, L. G. Wang, Y. Han, and Q. Y. He, “clusterProfiler: an R package for comparing biological themes among gene clusters,” OMICS, vol. 16, no. 5, pp. 284–287, 2012.
View at: Publisher Site | Google Scholar
D. P. Regan, L. Chow, S. Das et al., “Losartan blocks osteosarcoma-elicited monocyte recruitment, and combined with the kinase inhibitor toceranib, exerts significant clinical benefit in canine metastatic osteosarcoma,” Clinical Cancer Research, vol. 28, no. 4, pp. 662–676, 2022.
View at: Publisher Site | Google Scholar
Z. Hu, L. Li, W. Lan et al., “Enrichment of Wee1/CDC2 and NF-κB signaling pathway constituents mutually contributes to CDDP resistance in human osteosarcoma,” Cancer Research and Treatment, vol. 54, no. 1, pp. 277–293, 2022.
View at: Publisher Site | Google Scholar
M. Spalato and A. Italiano, “The safety of current pharmacotherapeutic strategies for osteosarcoma,” Expert Opinion on Drug Safety, vol. 20, no. 4, pp. 427–438, 2021.
View at: Publisher Site | Google Scholar
Y. Liu, Z. Zhang, Q. Li, L. Zhang, Y. Cheng, and Z. Zhong, “Mitochondrial APE1 promotes cisplatin resistance by downregulating ROS in osteosarcoma,” Oncology Reports, vol. 44, no. 2, pp. 499–508, 2020.
View at: Publisher Site | Google Scholar
T. Ando, J. Ichikawa, T. Fujimaki, N. Taniguchi, Y. Takayama, and H. Haro, “Gemcitabine and rapamycin exhibit additive effect against osteosarcoma by targeting autophagy and apoptosis,” Cancers, vol. 12, no. 11, p. 3097, 2020.
View at: Publisher Site | Google Scholar
C. W. Lee, M. C. Chi, T. M. Chang, and J. F. Liu, “Artocarpin induces cell apoptosis in human osteosarcoma cells through endoplasmic reticulum stress and reactive oxygen species,” Journal of Cellular Physiology, vol. 234, no. 8, pp. 13157–13168, 2019.
View at: Publisher Site | Google Scholar
X. Wang, X. Liang, H. Liang, and B. Wang, “SENP1/HIF-1α feedback loop modulates hypoxia-induced cell proliferation, invasion, and EMT in human osteosarcoma cells,” Journal of Cellular Biochemistry, vol. 119, no. 2, pp. 1819–1826, 2018.
View at: Publisher Site | Google Scholar
K. M. Goudarzi, M. Nistér, and M. S. Lindström, “mTOR inhibitors blunt the p53 response to nucleolar stress by regulating RPL11 and MDM2 levels,” Cancer Biology & Therapy, vol. 15, no. 11, pp. 1499–1514, 2014.
View at: Publisher Site | Google Scholar
S. Luo, J. Zhao, M. Fowdur, K. Wang, T. Jiang, and M. He, “Highly expressed ribosomal protein L34 indicates poor prognosis in osteosarcoma and its knockdown suppresses osteosarcoma proliferation probably through translational control,” Scientific Reports, vol. 6, no. 1, article 37690, 2016.
View at: Publisher Site | Google Scholar
D. D. Cheng, B. Zhu, S. J. Li, T. Yuan, Q. C. Yang, and C. Y. Fan, “Down-regulation of RPS9 inhibits osteosarcoma cell growth through inactivation of MAPK signaling pathway,” Journal of Cancer, vol. 8, no. 14, pp. 2720–2728, 2017.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Yang Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

251

Downloads

451

Citations