Next Article in Journal
Microvasculature Features Derived from Hybrid EPI MRI in Non-Enhancing Adult-Type Diffuse Glioma Subtypes
Next Article in Special Issue
Are CT-Derived Muscle Measurements Prognostic, Independent of Systemic Inflammation, in Good Performance Status Patients with Advanced Cancer?
Previous Article in Journal
EBV and Lymphomagenesis
Previous Article in Special Issue
Mucins as Potential Biomarkers for Early Detection of Cancer
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cancers
Volume 15
Issue 7
10.3390/cancers15072134
cancers-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Molecular Subtyping and Survival Analysis of Osteosarcoma Reveals Prognostic Biomarkers and Key Canonical Pathways

by
Siddesh Southekal
1,†,
Sushil Kumar Shakyawar
1,†,
Prachi Bajpai
2,
Amr Elkholy
2,
Upender Manne
2,
Nitish Kumar Mishra
1,* and
Chittibabu Guda
1,3,*
1
Department of Genetics, Cell Biology and Anatomy, University of Nebraska Medical Center, Omaha, NE 68198, USA
2
Department of Pathology, University of Alabama at Birmingham, Birmingham, AL 35233, USA
3
Center for Biomedical Informatics Research and Innovation, University of Nebraska Medical Center, Omaha, NE 68198, USA
*
Authors to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Cancers 2023, 15(7), 2134; https://0-doi-org.brum.beds.ac.uk/10.3390/cancers15072134
Submission received: 6 March 2023 / Revised: 28 March 2023 / Accepted: 30 March 2023 / Published: 4 April 2023
(This article belongs to the Collection Cancer Biomarkers)
Browse Figures
Review Reports Versions Notes

Abstract

:

Simple Summary

Osteosarcoma (OS), which accounts for about 35% of bone malignancy, shows aggressive progression in adults, adolescents, and children. Neoadjuvant chemotherapy provides an improvement in survival for OS patients; however, a molecular level understanding of the disease mechanisms is needed. We utilized publicly available multiomics data from the TARGET database, to perform survival analyses, differential expression analyses, mutational analyses, and subtyping using integrative clustering. Our results have identified several prognostic biomarkers (such as RAMP1, CRIP1, CORT, CHST13, and DDX60L) in OS patients that can be further explored for therapeutic applications.

Abstract

Osteosarcoma (OS) is a common bone malignancy in children and adolescents. Although histological subtyping followed by improved OS treatment regimens have helped achieve favorable outcomes, a lack of understanding of the molecular subtypes remains a challenge to characterize its genetic heterogeneity and subsequently to identify diagnostic and prognostic biomarkers for developing effective treatments. In the present study, global analysis of DNA methylation, and mRNA and miRNA gene expression in OS patient samples were correlated with their clinical characteristics. The mucin family of genes, MUC6, MUC12, and MUC4, were found to be highly mutated in the OS patients. Results revealed the enrichment of molecular pathways including Wnt signaling, Calcium signaling, and PI3K-Akt signaling in the OS tumors. Survival analyses showed that the expression levels of several genes such as RAMP1, CRIP1, CORT, CHST13, and DDX60L, miRNAs and lncRNAs were associated with survival of OS patients. Molecular subtyping using Cluster-Of-Clusters Analysis (COCA) for mRNA, lncRNA, and miRNA expression; DNA methylation; and mutation data from the TARGET dataset revealed two distinct molecular subtypes, each with a distinctive gene expression profile. Between the two subtypes, three upregulated genes, POP4, HEY1, CERKL, and seven downregulated genes, CEACAM1, ABLIM1, LTBP2, ISLR, LRRC32, PTPRF, and GPX3, associated with OS metastasis were found to be differentially regulated. Thus, the molecular subtyping results provide a strong basis for classification of OS patients that could be used to develop better prognostic treatment strategies.

Graphical Abstract

1. Introduction

Osteosarcoma (OS), a highly aggressive bone marrow malignancy that originates in mesenchymal tissue, accounts for about 35% of bone cancer cases worldwide. Other bone tumors include chondrosarcomas (25%), and Ewing sarcomas (16%) [1,2,3,4]. The most affected bones mainly include the humerus (10%), femur (43%), and tibia (23%), while sacrum areas, spine, and pelvis are rarely affected [5,6]. As per National Cancer Institute (NCI) data, for adults, OS (28%) is the second most prevalent primary bone cancer after chondrosarcoma (40%), while for children and adolescents, OS (56%) is the most common followed by Ewing sarcoma (34%). OS predominantly affects males [7,8]. Clinical manifestation suggests that, in 10–30% of the affected patients, metastases occur primarily in lungs (85%) and bones (8–10%) [5,9].
For various cancer types, including OS, genetic heterogeneity due to genomic instability poses challenges to the efficacy of current treatment strategies [10,11,12], demonstrating the need for accurate molecular subtyping to guide targeted therapeutic regimens for OS. Current advances in omics technologies, especially genomics, epigenomics, and transcriptomics, in combination with advanced data analytical tools have enabled molecular classification or subtyping of heterogeneous cancers that are increasingly preferred by physicians over the traditional and error-prone approaches based on clinical characteristics such as the histopathology, grade, and other visual observations [13]. To overcome the limitations of the previous approaches, recent methodologies have focused on the use of molecular profiling data for subtyping of cancers, including those of lung [14,15,16], colon [17,18], breast [19,20,21], and others [22,23,24]. In addition, new techniques such as NanoString and tissue microarray (TMA) approaches have also been evolved for subtype characterization of cancers [25,26,27,28].
In parallel, a variety of statistical and machine learning approaches such as Support Vector Machines (SVMs) and Deep Neural Networks (DNN) have been developed for molecular characterization and subclassification of cancers by assessing multiomics data [25,26]. These studies have been mainly fueled by the ease of access to public data repositories such as International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/, accessed on 1 December 2021), The Cancer Genome Atlas (TCGA), and Gene Expression Omnibus (GEO) for cancer sub-classification studies on the genomic, transcriptomic, and other clinical data.
In OS, the common subclassification method relies upon the characteristics of the dominant matrix present in the tumor. Conventionally, tumors occurring in fibrous tissue, bone, and cartilage are referred to as fibroblastic, osteoblastic, and chondroblastic OS, respectively [4,29]. These classes of tumors consist of about 75% of the cases; while the remaining 25% are considered to be variants, which show diverse biological properties. These variants are classified primarily based on the clinical factors, location of origin, and histological findings of the tumor. Similarly, the morphological and clinical feature-based classification includes parosteal, central low-grade, and periosteal OS. Tumors of these subclassifications constitute about 5% of the OS cases, but these have higher prognostic evaluations [30,31]. Despite an aggressive treatment regimen, the survival rate for OS is still a big challenge due to their propensity to metastasize before disease detection and to their resistance to therapies.
Unlike traditional histology-driven classifications, molecular subtyping of OS requires the identification of gene-based biomarkers and effective classifiers that can distinguish subtypes based on the molecular data. Prognostic biomarkers provide information about the patients’ overall cancer outcome and help to identify an appropriate therapeutic approach to cure or manage the disease. There is no systematic genome-wide analysis of these prognosis factors in OS [32,33,34,35]. Biomarker discovery in this area requires novel method development, which could improve the prognosis for OS patients. The current work focusses on developing novel gene-signature-based prognostic biomarkers in OS using multiomics data that include gene expression, DNA methylation, and miRNA expression data along with the associated clinical data. The findings may provide new insight in the domain of biomarker discoveries towards a better understanding of the subclassification of OS for developing improved treatment strategies.

2. Materials and Methods

Multiomic data: We downloaded gene expression, DNA methylation, miRNA expression, and clinical data for 85 common OS patients from Therapeutically Applicable Research To Generate Effective Treatments (TARGET, https://ocg.cancer.gov/, accessed on 1 December 2021) using a Bioconductor tool, TCGAbiolinks (an R/Bioconductor package for integrative analysis of TCGA data) [36]. We systematically accomplished data cleaning, global pattern analyses, and individual and integrative analyses on the multiomic datasets to remove noise and redundant data. Preprocessing steps of individual datasets are detailed in the following sections. Cox regression analysis and Kaplan Meier (KM) plots were generated to find genes associated with patient survival. Patient samples common across all data types were considered for downstream clustering analysis. Pathway analysis was carried out to identify the functional relevance of the differentially expressed genes between the identified molecular subtypes in OS.
Protein coding and lncRNA data: TCGAbiolinks was used to query and download the HTSeq–FPKM (Fragments Per Kilobase of Exon per Million) gene expression count data from the TARGET-OS project. Genes having more than 20% missing values across patients were not considered further for analysis. Additionally, only those genes with CPM (count per million) values > 1 for more than 10 % samples were retained. Genes were annotated as protein-coding or lncRNA using gencode reference annotation v22. The normalized count data were gene median centered, and the top 2500 protein-coding and 500 lncRNA based on median absolute deviation (MAD) were retained for further clustering analysis.
miRNA data: Normalized miRNA data were downloaded from NCI TARGET’s web portal. The dataset was first filtered by removing duplicated miRNA expression data. Further, clustering was performed using nonnegative matrix factorization (NMF) by retaining 190 (25%) of the most variable miRNAs.
Methylation array data: The raw IDAT files from the methylation study (NCBI GEO accession: GSE72872, Illumina HumanMethylation450 BeadChip array) [37] were downloaded using the R function, getGEOSuppFiles. The IDAT files were processed using the openSesame [38] tool in R to obtain the intensities of methylated and unmethylated alleles at the CpG sites as the beta value numeric. CpG probes with missing values in more than 20% samples across the patients were removed from further analysis. Imputation was performed to fill the missing β values in the data matrix using the imputeKNN function, with number of neighbor (k) of 15 [39]. This was followed by beta mixed-integer quantile (BMIQ) normalization using the ChAMP tool [40]. CpG probes mapped against sex and mitochondrial chromosomes were excluded from analyses to eliminate gender bias. CpG probes that overlapped with repeat masker and SNPs from dbSNP v151 with minor allele frequency (MAF) greater than 1% were removed to draw out sequence polymorphisms. For the final matrix, 5000 probes were retained based on the standard deviation for further downstream clustering analysis.
Mutation data: Mutation data from 123 OS patients was downloaded as a tab-delimited Mutation Annotation Format (MAF) file along with the scores, all-lesions, amplification, and deletion gene files generated by the Genepattern GISTIC module [35]. These data were used as input for the Bioconductor package, maftools [36]. Genes having a mutation in less than 5% of the samples were removed from further analysis. The matrix was converted into binary form for the presence or absence of mutations. The top 25% most variable genes (n = 1107) having a mutation in more than 20% of the samples were retained in the final matrix for clustering analysis.
Clustering analysis: NMF based clustering was performed using nrun = 500 for the number of clusters, k = 2 to 10 for miRNA, mutation, protein-coding, lncRNA, and methylation data for the same set of patients across all datasets by using R package NMF. The preferred clustering result was determined by using the observed cophenetic correlation coefficients between clusters and the average silhouette width of consensus cluster members.
Cluster-Of-Clusters Analysis (COCA): The R package, COCA [41], an integrative clustering method was used to summarize clusters found in the protein-coding, lncRNA, miRNA, methylation, and mutation datasets by identifying a global clustering that is in good agreement with the clustering output in each of the individual datasets. The final clustering was accomplished based on the construction of the Matrix-Of-Clusters (MOC), a binary matrix of size N × K, where K is the sum of the number of clusters.
Differential gene expression analysis: Differential gene expression (DGE) analysis between the samples belonging to the two molecular subtypes obtained from the COCA was performed using the Bioconductor tool, DESeq2. The Benjamini-Hochberg (BH) adjusted p-value cut-off of 0.01, and an absolute log2Fold Change of 1 was used to obtain the list of differentially expressed genes. The heatmap was generated using tidyheatmap package in R.
Pathway analysis: All differentially expressed genes (with adjusted p-value < 0.01 and log2Fold Change = 1) between Subtypes 1 and 2 were used as input into QIAGEN IPA (QIAGEN, Redwood City, CA, USA; www.qiagen.com/ingenuity, accessed on 1 December 2021) for pathway enrichment analysis.
Survival analysis: We performed survival analyses for protein-coding, lncRNA, and miRNA genes using the log-rank test and Cox-regression analysis (p-value ≤ 0.05) by constructing high and low expression groups, using the median expression of genes as the cut-off value. For gene expression analysis, we used HTSeq–FPKM gene expression count data of OS patients. The R tool, “survival” was used for survival analysis, and the Kaplan–Meier survival curve plots were generated for all analyses. In addition, a log-rank or Mantel–Haenszel test was conducted to assess differences between survival curves of the two groups and to calculate the p-values. A hazard ratio (HR) > 1 indicates that patients in the high expression group had low survival, and HR < 1 suggested high or better survival.
Correlation analysis: Pearson correlation analysis between gene expression for the protein-coding genes was performed using the R cor.test function. The association was considered significant at a p-value < 0.05.

3. Results

In the present work, we have utilized the OS omics data of five distinct datatypes which include miRNA, mRNA, lncRNA, methylation, and mutation data, along with the survival data from public resources including the TARGET database and the PubMed repository. Initially, three data types (miRNA, mRNA, and lncRNA) along with the survival information were used to understand the associations of these biomolecules to the survival of OS patients. The later part of the study mainly focused on COCA-based clustering of 85 patients using multi-level omics information including miRNA, mRNA, lncRNA, methylation, and mutation, resulting in the identification of two OS subtypes. Results from each analysis are described in the following sections.

3.1. Mucin Family Genes Are Highly Mutated in OS

For mutational signature analysis, we read MAF files using maftool, which provided frequencies of six classes of the mutation types (C > A, C > G, C > T, T > A, T > C, and T > G) in each sample. The contribution of each mutational signature was analyzed across the samples. Our analyses showed that the commonly mutated genes include PABPC3, HRNR, KMT2C, and PRSS3 in more than 65% of the samples (Supplementary Figure S1A). A missense mutation was found to be the most common variant type across all patients, although many of the genes were found to have multiple mutations as Multi_Hit (green), as indicated in Supplementary Figure S1A. The analysis also showed that the top-ranked altered genes include the MUC family of genes, namely, MUC6 (80%), MUC12 (76%), and MUC4 (76%). An oncoplot of members of the MUC family genes with their mutation frequencies is shown in Supplementary Figure S1B.
Enriched pathways associated with OS genetic alterations were analysed using the OncogenicPathways module in maftools for known signaling pathways in TCGA cancers [42]. These include the P53, Notch, Myc, Wnt, Hippo, RTK-RAS, PI3K, and Cell cycle signaling pathways. The top five most commonly altered genes prevalent in more than 10% of the samples containing mutations in the TCGA enriched pathways are shown in Figure 1A. The genes associated with these pathways include TEAD2, WWC1, AJUBA, CSNK1D, and LLGL2 in hippo signaling, IGF1R, RAC1, ERBB4, GRB2, and NF1 in RTK-RAS pathway, RFNG, LFNG, MAML3, APH1A, and NOTCH2 in Notch pathway, and WNT7B, LZTR1, DKK1, APC, and FZD6 in Wnt pathway. In the TP53 pathway, genes including CHEK2, and ATM mostly showed missense mutations (blue color), while the TP53 gene had multiple mutations including missense, Amplification (Amp), Deletion (Del), and Frameshift. NOTCH signaling pathways mostly showed Amp mutation as presented in brick color. Similarly, the genes LLGL2, CSNK1D, and AJUBA in the HIPPO pathway showed Amp mutation; while other genes, WWC1 and TEAD2 in the same pathway, were deleted. Cell cycle pathways include genes RB1 which showed multiple mutation types, whereas, other genes CCNE1 and CDK4 showed mainly amplification. Multiple types of mutation were identified in RB1 in Cell Cycle and TP53 in TP53 pathways, NAML3 in NOTCH pathways, MGA in MYC pathways, and NF1 in RTK-RAS pathways. In general, amplifications were most prevalent across pathways. The fraction of affected genes in the most enriched pathways included 94% in Hippo, 86% in Wnt, 84% in RTK-RAS, and 84% in Notch signaling pathways (Figure 1B). The pathways, associated genes, fraction of mutated samples, and other related information is provided in Supplementary Table S1.

3.2. Survival Analysis Reveals Key Prognostic Pathway Genes

Survival analysis of OS patients was performed using gene expression, miRNA, and lncRNA, with p-values < 0.05 for both log-odd and Cox regression. In-house R code with survival and survminer packages in the background were used to perform this analysis.
The log-Rank test and cox-regression analysis (p-value < 0.05) using protein-coding HTSeq-FPKM gene expression identified 430 genes associated with OS prognosis (p-value < 0.01). The genes with HR more than 4.5 include RAMP1, TAC4, MT1A, CRIP1, NHEJ1, TIMM23B, COL5A2, UFC1, and CORT. The higher expression of these genes indicate low survival of OS patients. In contrast, higher expression of genes such as METTL20, MYH10, PDE1B, CCDC42, and many more which got HR < 1 refers to better survival rates. In each category, survival plots for the top three potential genes are provided in Figure 2A,B, respectively. The HR and p-value of all 430 genes are provided in Supplementary Table S2A subset of the top 49 genes, significantly (p-value < 0.001) related to the survival of OS patients, are shown in Table 1.
We mapped the expression of 34 genes involved in 10 canonical pathways that are considered to be the likely cancer drivers (functional contributors) or therapeutic targets [42], and are found to be associated with the survival of patients (p < 0.05) in our analysis. As shown in Figure 3A, most of the genes had either positive or no correlation within the pathways at p < 0.05. Only a few genes such as SOST and WNT5A in WNT signaling pathways, MAPK1 and PEBP1 in RTK-RAS pathways, and PIK3R3 and RPS6 in PI3K pathways, showed negative expression pattern within the pathways. In our analyses, none of the genes in the NRF2 and TGFβ signaling pathways were significant in survival and are, therefore, not shown in Figure 3A. HR values of 34 genes from these pathways are shown in Figure 3B. Genes including ATM in TP53 pathway, CDKN2A in Cell cycle, FAT1 and FAT3 in HIPPO signaling, MXl1 in MYC pathway, DTX3, DTX2, JAG2, HES4, and HES5 in NOTCH signaling showed HR values greater than 1 indicating their association with low survival in OS patients. Other genes such as LRP5, SOST, WNT1, RPS6 in different pathways also have HR value more than 1 with possible association with survival (Figure 3B). In contrast, genes such as KRAS, WNT5A, and CDK6 have a HR less than 1 indicating their association with better survival in OS patients.
Survival analysis in relation to lncRNAs showed that 36 lncRNAs (p-value < 0.01) or 142 lncRNAs (p-value < 0.05) are correlated with survival. ELFN1-AS1 (HR > 1) is among the top 3 lncRNAs identified (p-value < 0.001) in our study, which has been recognized as an important signature to predicting OS patient survival in a previous study [43]. Similarly, RP11-283C24.1 and lnc01060 also displayed HR > 1, indicating that higher expression of these lncRNAs may lead to lower survival of OS patients. No information was found in the literature about two other lncRNAs, lnc01060 and RP11-283C24.1 and further investigations are required to validate these potential biomarkers. Our analysis also identified lncRNAs CTD-2078B5.2, RP11-78O7.2, and RP11-446N19.1, which got HR < 1 and p-value < 0.05. Higher expression of these lncRNAs are associated with better survival of the patients; hence, they serve as potential biomarkers for OS prognosis. Survival plots of the top lncRNAs with HR > 1 and HR < 1 are provided in the Supplementary Figure S3. The HR and beta values of all 36 and 142 lncRNAs identified above are provided in the Supplementary Table S3.
In addition, survival analysis using miRNA data resulted in the identification of 76 (p < 0.05) or 31 miRNAs (p < 0.01), which showed association with the survival (Supplementary Table S4). We identified four miRNAs, miR-122, miR-200b, miR-1298, and miR-1264 with HR < 1 and p-value < 0.05. Similarly, miRNAs with HR > 1 and p-value < 0.05 include hsa-miR-182, hsa-miR-891a, hsa-miR-190, miR-452, hsa-miR-142, and 62 more. Similar to mRNAs and lncRNAs, higher expression of miRNAs with HR > 1 indicate lower survival of the patients; whereas, higher expression of miRNAs with HR < 1 shows better survival association with the OS patients. Previous studies have shown that miRNA, hsa-miR-190 was involved in many cancer-related biological processes such as metastasis, proliferation, and apoptosis by means of dysregulating target genes [44,45,46]. Similarly, hsa-miR-452 and hsa-miR-122 were found to be involved in promoting colorectal and colon cancer progression, respectively [47,48]. The Kaplan-Meier survival plots for top miRNAs are shown in Supplementary Figure S4.

3.3. Clustering Analysis

Clustering analysis with NMF was performed for 85 patients using the miRNA, mutation, protein-coding, lncRNA, and CpG probe methylation datasets. The final matrix consisted of 190 miRNAs, 1107 mutations, 2500 protein-coding regions, 500 lncRNA, and 5000 CpG methylation sites, which were obtained after the preprocessing steps as described previously in the methods section. The optimal number of clusters was obtained by evaluating the cophenetic correlation coefficients between clusters and the average silhouette width (ASW), resulting in three clusters each for protein-coding, lncRNA, miRNA, and methylation data, and four clusters for the mutation data (Figure 4).
A matrix of clusters was generated and used as input for COCA, which resulted in two major integrated molecular clusters/subtypes of OS, consisting of 44 and 41 patients in Subtype-1 and Subtype-2, respectively. The matrix of clusters was plotted using plotMOC function. Supplementary Figure S2.

3.4. Comparison of the Identified Molecular Subtypes

Differential Gene Expression analysis was performed between the two molecular subtypes obtained from COCA using the Bioconductor package, DESeq2 [49]. A total of 261 genes were differentially expressed between Subtype-2 vs. Subtype-1 at log2Fold Change > 1.5 and adjusted p-value < 0.05 (Figure 5A). Among these, 178 were upregulated and 83 were downregulated. A heatmap of the top 100 most variables differentially expressed genes is shown in Figure 5B. Hierarchical clustering based on gene expression data showed two distinguishable molecular subtypes consistent with the multiomics-based subclusters. The expression-based heatmap (Figure 5B) clearly showed the differences in gene expression levels between Subtype 1 (blue color) and Subtype 2 (red color). We utilized gene count data from DEG’s for clustering using euclidian distance measure for hierarchical clustering with tidyheatmap library in R, which also resulted in two major clusters consistent with the multiomics-based molecular subtyping (Supplementary Table S5, Figure 5B). However, no genes were found to be significantly differentially mutated between the two molecular subtypes. The pathway enrichment analysis of survival associated protein-coding genes (as identified in log-odd test and Cox-regression analysis, p-value < 0.05) using IPA identified WNT signaling and Axonal Guidance signaling pathways, which are related to the OS pathogenesis [50,51,52] (Supplementary Figure S5). However, no significant difference was observed in the survival of patients between the two subtypes.

4. Discussion

OS is one of the most frequent and aggressive malignant bone neoplasms. Identification of clinically relevant prognostic markers in OS is of vital importance. Cancer subtyping studies using multiomics data have been widely performed to get an understanding of the genetic mechanism of the disease as well as to accelerate their therapeutic applications [53,54,55,56]. Although, many markers have been demonstrated to be of prognostic significance in the treatment of OS, there is a need for new therapeutic targets that can be evaluated [33,57]. In the current study, we accomplished global analysis of different types of omics data from the OS patients in the TARGET database, and identified molecular subtypes by integrative analysis of multiomics data.
Mutation analysis identified significant mutations in the Mucin (MUC) family of genes, notably MUC6 (80% of the patients), MUC12 (76% of the patients), and MUC4 (76% of the patients). Previously, Tirabosco et al. have reported MUC4 as the immunoreactive gene in fibromyxoid and sclerosing tumors in bone [58]. Similarly, MUC12 and MUC4 were identified as top-altered genes in a cohort of 21 OS patients [59]. However, the understanding of whether or not these mutations in MUC genes are driving OS is still unknown. In our analyses, MUC15 which promote cell proliferation, migration, and invasion [60], mutated in only 2% of the TARGET OS patient cohort.
Further, the gene KMT2C, which is mutated in 67% of all OS samples, involved in OS carcinogenesis and progression [61]. Similar to results of our study, the exome analysis shows a high load of somatic variations in KMT2C, which can have histone modification-related activity in OS [62]. Further investigation about the roles of identified mutations in various genes could facilitate the development of prognostic biomarkers for OS. Subsequent pathway-level analysis showed that the Wnt pathway is one of the top enriched pathways in OS samples and is significantly involved in OS progression, as supported by several previous studies [50,63,64].
We also observed enrichment of pathway members that contained recurrent or known alterations that are likely to be cancer drivers or therapeutic targets in the Hippo, RTK-RAS, Notch, and Wnt signaling pathways [42] (Figure 1). The Hippo and Notch signaling pathways have an active role in OS [65,66,67]. Fang et al., 2018 have shown that the hyperactive Wnt/β-catenin pathway may be required for proliferation and metastasis of OS [68]. In addition, similar to the present results, somatic mutations and upregulated gene expression of many components in the Wnt/β-catenin pathway were observed in their study.
Survival analysis revealed the survival association of several genes to OS including those with high expression such as the receptor activity-modifying proteins 1 (RAMP1): a protein of calcium signaling receptor complex; Tachykinin Precursor 4 (TAC4), Metallothionein 1A (MT1A), CRIP1, and CORT. Our analysis showed that these genes are associated with low survival of OS patients based on their HR values (RAMP1: 5.7, TAC4: 5.5, MT1A: 4.8, CRIP: 4.8, CORT: 4.6) (Figure 2B, Supplementary Table S2). Similarly, CHST13 (HR: 2.4) is associated with low survival, while DDX60L (HR 0.37), METTL20 (HR 0.14), MYH10 (HR 0.16), and PDE1B (HR 0.15) are associated with higher survival outcomes in OS patients (Figure 2A, Supplementary Table S2). These genes have the potential to serve as prognostic biomarkers in OS. The prognostic role of RAMP1 in OS has not been reported in the literature. RAMP1 interacts with G-protein-coupled receptors, such as calcitonin receptor-like receptor, calcitonin receptor, calcium-sensing receptor, and glucagon receptor [69,70]. Knockdown of RAMP1 reduced clonogenic/spheroidal growth and tumorigenicity, and small-molecule inhibitors directed against the RAMP1 reduced growth of Ewing sarcoma [69].
A lower survival rate in patients was related to high TAC4 and MT1A expression. Although the role of TAC4 in metastasis and prognosis of OS is unknown, other bioinformatic studies have shown a possible link [71]. MT1A expression levels were used as one of the genes to develop a machine learning-based prognostic risk model in OS patients, suggesting that this gene is associated with OS survival [72].
From our study, higher expression of CRIP1 and CORT was associated with lower survival (p < 0.001). CRIP1 has been previously shown to be over-expressed in pre-therapeutic OS samples [73]. CORT is an endogenous cyclic neuropeptide known to regulate the growth and metastasis of lung and thyroid cancer [74,75]. The higher expression of CORT gene has also been associated with high-risk group OS patients [76]. In addition, we have also identified lncRNA and miRNA which were associated with OS survival. In particular, higher expression of miRNAs hsa-miR-190, hsa-miR-452, and hsa-miR-142 are directly associated with lower survival of the OS patients. In contrast, hsa-miR-122, hsa-miR-200, and hsa-miR-1298 had association with improved survival of the patients. Similarly, lncRNAs ELFN1-AS1, lnc01060, RP11-283C24.1, and RP11-283C24.1 with HR > 1 and lncRNAs CTD-2078B5.2, RP11-78O7.2, and RP11-446N19.1 with HR < 1 were interpreted to be associated with lower and higher survival of the OS patients. The experimental validation is required to support these findings.
Our global clustering strategy built upon the multiomics landscape has identified two main subtypes (Subtype 1: 44 patients, Subtype 2: 41 patients), which opens up the potential for developing effective treatments. Differential expression and pathway enrichment analyses suggest that both subtypes had different deregulated pathways associated with OS pathology. We observed over 250 differentially expressed genes between two subtypes (log2Fold Change > 1.5 and adjusted p-value < 0.05), including genes associated with Wnt signaling, calcium signaling, and PI3K-Akt signaling pathways, were highly enriched and found related to OS pathogenesis.
Although no difference in somatic mutation and somatic copy number alteration levels between the identified subtypes, hierarchical clustering of the top 100 most differentially expressed genes between the two OS subtypes could be used to classify the patients into two molecular clusters. Liu et al. 2020 [1] identified differentially regulated genes i.e., VAMP8, A2M, HLA-DRA, SPARCL1, HLA-DQA1, APOC1, and AQP1, which are involved in metastasis of OS. Our analysis showed that genes VAMP8, SPARCL1, and APOC1 were significant with p-value 0.007, 0.017, and 0.038, respectively, for predicting survival in OS patients. Our analysis also identified three upregulated genes, POP4, HEY1, CERKL, and seven other downregulated genes, CEACAM1, ABLIM1, LTBP2, ISLR, LRRC32, PTPRF, and GPX3 between identified two subtypes.

5. Conclusions

This study has identified potential prognostic genes including CRIP1, CORT, CHST13, and DDX60L in the OS samples. Although many other cancer-related genes showed a high frequency of mutations in OS, the MUC family of genes (MUC6, MUC12, and MUC4) were found to be highly mutated in our study. miRNAs: hsa-miR-190, hsa-miR-452, hsa-miR-142, hsa-miR-122, and hsa-miR-200, and lncRNAs: ELFN1-AS1, lnc01060, RP11-283C24.1, RP11-283C24.1, CTD-2078B5.2, RP11-78O7.2, and RP11-446N19.1 were also found to be significantly associated with the survival of OS patients; and therefore, interpreted as potential biomarkers that can to be further validated with experimental studies. Our detailed clustering approach using COCA has identified two main OS subtypes, which are significantly different at the molecular level. These findings offer a strong basis to understand the genetic heterogeneity of OS and develop next generation prognostic biomarkers for effective treatment of OS patients.

Supplementary Materials

The following supporting information can be downloaded at: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/cancers15072134/s1, Figure S1: (A) Oncoplot of the top 15 most frequently mutated genes 123 OS patients. Colored squares show mutated genes, while grey squares show non mutated genes. Each color represents a different type of mutations. Each color represents different type of mutations as labelled in the figure. (B) Oncoplot of the frequency of mutation in MUC family genes in 123 OS patients. Colored squares show mutated genes, while grey squares show non mutated genes. Each color represents a different type of mutations. Each color represents different type of mutations as labelled in figure; Figure S2: Clustering of 85 OS patients into two molecular subtypes based on multiomics data using cluster-of-cluster analysis. The final multiomics-based molecular subtypes consisted of 44 and 41 samples respectively. The vertical sidebar shows the annotation of cluster labels obtained using individual of ‘omics’ data and the final multiomics-based clustering for each sample; Figure S3: Kaplan-Meier plot for high versus low expression group in TARGET OS data with p-value from log-rank test and Cox regression model for lncRNAs CTD-2078B5.2, RP11-78O7.2, and RP11-446N19.1 (HR < 1) (A), and ELFN1-AS1, RP11-283C24.1 and lnc01060 (HR > 1) (B); Figure S4: Kaplan-Meier plot for high vs low expression group in TARGET OS data with p-value from log-rank test and Cox regression model for miRNAs miR-122, miR-200b, miR-1298 (HR < 1) (A), and miR-190, miR-452, and miR-488 (HR > 1) (B); Figure S5: Top canonical pathway enriched between OS subtype 1 and subtype 2. (Analysis performed by IPA). Table S1: Pathways, associated genes and fraction of mutated samples in each pathway; Table S2: Kaplan-Meier plot for high vs low expression group in TARGET OS data with p-value from log-rank test and Cox regression model for genes METTL20, MYH10, and PDE1B with HR < 1 (A), and genes RAMP1, TAC4, and MT1 with HR > 1(B); Table S3: List of lncRNA which are significantly associated with the survival of OS patients. The Hazard Ratio of each lncRNA along with p-values and other information is also mentioned; Table S4: List of miRNAs which are significantly associated with the survival of OS patients. The Hazard Ratio of each miRNA along with p-values and other information is also mentioned; Table S5: Results from differential gene expression analysis between Subtype 1 and Subtype 2. Genes with corresponding fold change and adjusted p-value, and other related information are listed.

Author Contributions

Conceptualized and designed the work plan, S.S., N.K.M. and C.G.; performed the analysis and produced the figures, S.S.; performed genes and pathways enrichment, and survival analysis, S.K.S.; wrote and revised the manuscript, S.S. and S.K.S.; assisted in survival analysis and with manuscript writing, P.B., A.E. and U.M.; edited and improved the manuscript, N.K.M. and C.G. All authors have read and agreed to the published version of the manuscript.

Funding

We acknowledge the support from Nebraska Research Initiative (NRI) and NIH awards (2P01AG029531, 2P20GM103427, 5P30CA036727, 2U54GM115458) to CG, and the grant (P30CA013148) funded to the UAB Tissue Biorepository Facility of the UAB O’Neal Comprehensive Cancer Center.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All analyses were performed using the R version 4.0.1 (R Development Core Team 2015). We performed clustering, pathway and survival analysis by using R/-Bioconductor tools. All the relevant data generated in the project is made available in the supplementary files.

Acknowledgments

The authors would like to express their gratitude to the Bioinformatics and Systems Biology Core at UNMC for providing the computational infrastructure, the Holland Computing Center (HCC) for offering access to supercomputers. The authors also thank Neetha Vellichirammal for providing valuable suggestions on data interpretation and revising manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lancia, C.; Anninga, J.K.; Sydes, M.R.; Spitoni, C.; Whelan, J.; Hogendoorn, P.C.W.; Gelderblom, H.; Fiocco, M. A novel method to address the association between received dose intensity and survival outcome: Benefits of approaching treatment intensification at a more individualised level in a trial of the European Osteosarcoma Intergroup. Cancer Chemother. Pharmacol. 2019, 83, 951–962. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Menéndez, S.; Gallego, B.; Murillo, D.; Rodríguez, A.; Rodríguez, R. Cancer Stem Cells as a Source of Drug Resistance in Bone Sarcomas. J. Clin. Med. 2021, 10, 2621. [Google Scholar] [CrossRef] [PubMed]
  3. Cortini, M.; Baldini, N.; Avnet, S. New Advances in the Study of Bone Tumors: A Lesson From the 3D Environment. Front. Physiol. 2019, 10, 814. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Misaghi, A.; Goldin, A.; Awad, M.; Kulidjian, A.A. Osteosarcoma: A comprehensive review. SICOT-J 2018, 4, 12. [Google Scholar] [CrossRef] [Green Version]
  5. Bielack, S.S.; Kempf-Bielack, B.; Delling, G.; Exner, G.U.; Flege, S.; Helmke, K.; Kotz, R.; Salzer-Kuntschik, M.; Werner, M.; Winkelmann, W.; et al. Prognostic Factors in High-Grade Osteosarcoma of the Extremities or Trunk: An Analysis of 1,702 Patients Treated on Neoadjuvant Cooperative Osteosarcoma Study Group Protocols. J. Clin. Oncol. 2002, 20, 776–790. [Google Scholar] [CrossRef] [PubMed]
  6. Luetke, A.; Meyers, P.A.; Lewis, I.; Juergens, H. Osteosarcoma treatment—Where do we stand? A state of the art review. Cancer Treat. Rev. 2014, 40, 523–532. [Google Scholar] [CrossRef]
  7. Jo, V.Y.; Fletcher, C.D. WHO classification of soft tissue tumours: An update based on the 2013 (4th) edition. Pathology 2014, 46, 95–104. [Google Scholar] [CrossRef]
  8. Sbaraglia, M.; Bellan, E.; Tos, A.P.D. The 2020 WHO Classification of Soft Tissue Tumours: News and perspectives. Pathologica 2020, 113, 70–84. [Google Scholar] [CrossRef]
  9. Odri, G.A.; Tchicaya-Bouanga, J.; Yoon, D.J.Y.; Modrowski, D. Metastatic Progression of Osteosarcomas: A Review of Current Knowledge of Environmental versus Oncogenic Drivers. Cancers 2022, 14, 360. [Google Scholar] [CrossRef]
  10. Burrell, R.A.; McGranahan, N.; Bartek, J.; Swanton, C. The causes and consequences of genetic heterogeneity in cancer evolution. Nature 2013, 501, 338–345. [Google Scholar] [CrossRef]
  11. Wang, D.; Niu, X.; Wang, Z.; Song, C.-L.; Huang, Z.; Chen, K.-N.; Duan, J.; Bai, H.; Xu, J.; Zhao, J.; et al. Multiregion Sequencing Reveals the Genetic Heterogeneity and Evolutionary History of Osteosarcoma and Matched Pulmonary Metastases. Cancer Res. 2019, 79, 7–20. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Liu, W.; Wang, R.; Zhang, Y.; Wang, H.; Huang, Z.; Jin, T.; Yang, Y.; Sun, Y.; Cao, S.; Niu, X. Whole-exome sequencing in osteosarcoma with distinct prognosis reveals disparate genetic heterogeneity. Cancer Genet. 2021, 256, 149–157. [Google Scholar] [CrossRef]
  13. Zheng, Y.; Huang, Y.; Bi, G.; Du, Y.; Liang, J.; Zhao, M.; Chen, Z.; Zhan, C.; Xi, J.; Wang, Q. Multi-omics characterization and validation of MSI-related molecular features across multiple malignancies. Life Sci. 2021, 270, 119081. [Google Scholar] [CrossRef] [PubMed]
  14. Niemira, M.; Collin, F.; Szalkowska, A.; Bielska, A.; Chwialkowska, K.; Reszec, J.; Niklinski, J.; Kwasniewski, M.; Kretowski, A. Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers 2019, 12, 37. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Schwendenwein, A.; Megyesfalvi, Z.; Barany, N.; Valko, Z.; Bugyik, E.; Lang, C.; Ferencz, B.; Paku, S.; Lantos, A.; Fillinger, J.; et al. Molecular profiles of small cell lung cancer subtypes: Therapeutic implications. Mol. Ther. Oncolytics 2021, 20, 470–483. [Google Scholar] [CrossRef] [PubMed]
  16. Wang, X.; Yu, G.; Wang, J.; Zain, A.M.; Guo, W. Lung cancer subtype diagnosis using weakly-paired multi-omics data. Bioinformatics 2022, 38, 5092–5099. [Google Scholar] [CrossRef]
  17. Valenzuela, G.; Canepa, J.; Simonetti, C.; de Zaldívar, L.S.; Marcelain, K.; González-Montero, J. Consensus molecular subtypes of colorectal cancer in clinical practice: A translational approach. World J. Clin. Oncol. 2021, 12, 1000–1008. [Google Scholar] [CrossRef]
  18. Zhou, B.; Yu, J.; Cai, X.; Wu, S. Constructing a molecular subtype model of colon cancer using machine learning. Front. Pharmacol. 2022, 13, 1008207. [Google Scholar] [CrossRef]
  19. Eliyatkin, N.; Yalcin, E.; Zengel, B.; Aktaş, S.; Vardar, E. Molecular Classification of Breast Carcinoma: From Traditional, Old-Fashioned Way to A New Age, and A New Way. J. Breast Health 2015, 11, 59–66. [Google Scholar] [CrossRef] [Green Version]
  20. Fragomeni, S.M.; Sciallis, A.; Jeruss, J.S. Molecular Subtypes and Local-Regional Control of Breast Cancer. Surg. Oncol. Clin. N. Am. 2018, 27, 95–120. [Google Scholar] [CrossRef]
  21. Ochoa, S.; Hernández-Lemus, E. Functional impact of multi-omic interactions in breast cancer subtypes. Front. Genet. 2023, 13, 1078609. [Google Scholar] [CrossRef]
  22. Martini, G.; Dienstmann, R.; Ros, J.; Baraibar, I.; Cuadra-Urteaga, J.L.; Salva, F.; Ciardiello, D.; Mulet, N.; Argiles, G.; Tabernero, J.; et al. Molecular subtypes and the evolution of treatment management in metastatic colorectal cancer. Ther. Adv. Med Oncol. 2020, 12, 1758835920936089. [Google Scholar] [CrossRef] [PubMed]
  23. Alizadeh, A.A.; Eisen, M.B.; Davis, R.E.; Ma, C.; Lossos, I.S.; Rosenwald, A.; Boldrick, J.C.; Sabet, H.; Tran, T.; Yu, X.; et al. Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature 2000, 403, 503–511. [Google Scholar] [CrossRef] [PubMed]
  24. Collisson, E.A.; Sadanandam, A.; Olson, P.; Gibb, W.J.; Truitt, M.; Gu, S.; Cooc, J.; Weinkle, J.; Kim, G.E.; Jakkula, L.; et al. Subtypes of pancreatic ductal adenocarcinoma and their differing responses to therapy. Nat. Med. 2011, 17, 500–503. [Google Scholar] [CrossRef] [PubMed]
  25. Inaki, K.; Shibutani, T.; Maeda, N.; Eppenberger-Castori, S.; Nicolet, S.; Kaneda, Y.; Koyama, K.; Qiu, Y.; Wakita, K.; Murakami, M. Pan-cancer gene expression analysis of tissue microarray using EdgeSeq oncology biomarker panel and a cross-comparison with HER2 and HER3 immunohistochemical analysis. PLoS ONE 2022, 17, e0274140. [Google Scholar] [CrossRef] [PubMed]
  26. He, J.; Chen, Z.; Xue, Q.; Sun, P.; Wang, Y.; Zhu, C.; Shi, W. Identification of molecular subtypes and a novel prognostic model of diffuse large B-cell lymphoma based on a metabolism-associated gene signature. J. Transl. Med. 2022, 20, 186. [Google Scholar] [CrossRef] [PubMed]
  27. Chen, Y.-J.; Huang, C.-S.; Phan, N.-N.; Lu, T.-P.; Liu, C.-Y.; Huang, C.-J.; Chiu, J.-H.; Tseng, L.-M. Molecular subtyping of breast cancer intrinsic taxonomy with oligonucleotide microarray and NanoString nCounter. Biosci. Rep. 2021, 41. [Google Scholar] [CrossRef]
  28. Kulkarni, M.M. Digital Multiplexed Gene Expression Analysis Using the NanoString nCounter System. Curr. Protoc. Mol. Biol. 2011, 94, 25B-10. [Google Scholar] [CrossRef]
  29. Raymond, A.K.; Jaffe, N. Osteosarcoma Multidisciplinary Approach to the Management from the Pathologist’s Perspective. Pediatr. Adolesc. Osteosarcoma 2009, 152, 63–84. [Google Scholar] [CrossRef]
  30. Zhu, L.; McManus, M.M.; Hughes, D.P.M. Understanding the Biology of Bone Sarcoma from Early Initiating Events through Late Events in Metastasis and Disease Progression. Front. Oncol. 2013, 3, 230. [Google Scholar] [CrossRef] [Green Version]
  31. Kurt, A.M.; Unni, K.K.; McLeod, R.A.; Pritchard, D.J. Low-grade intraosseous osteosarcoma. Cancer 1990, 65, 1418–1428. [Google Scholar] [CrossRef] [PubMed]
  32. Zhao, N.; Guo, M.; Wang, K.; Zhang, C.; Liu, X. Identification of Pan-Cancer Prognostic Biomarkers Through Integration of Multi-Omics Data. Front. Bioeng. Biotechnol. 2020, 8, 268. [Google Scholar] [CrossRef]
  33. Zamborsky, R.; Kokavec, M.; Harsanyi, S.; Danisovic, L. Identification of Prognostic and Predictive Osteosarcoma Biomarkers. Med. Sci. 2019, 7, 28. [Google Scholar] [CrossRef] [Green Version]
  34. Shi, Y.; He, R.; Zhuang, Z.; Ren, J.; Wang, Z.; Liu, Y.; Wu, J.; Jiang, S.; Wang, K. A risk signature-based on metastasis-associated genes to predict survival of patients with osteosarcoma. J. Cell. Biochem. 2020, 121, 3479–3490. [Google Scholar] [CrossRef]
  35. Qi, W.; Yan, Q.; Lv, M.; Song, D.; Wang, X.; Tian, K. Prognostic Signature of Osteosarcoma Based on 14 Autophagy-Related Genes. Pathol. Oncol. Res. 2021, 27, 1609782. [Google Scholar] [CrossRef] [PubMed]
  36. Colaprico, A.; Silva, T.C.; Olsen, C.; Garofano, L.; Cava, C.; Garolini, D.; Sabedot, T.S.; Malta, T.M.; Pagnotta, S.M.; Castiglioni, I.; et al. TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016, 44, e71. [Google Scholar] [CrossRef] [PubMed]
  37. Krause, L.; Nones, K.; Loffler, K.A.; Nancarrow, D.; Oey, H.; Tang, Y.H.; Wayte, N.J.; Patch, A.M.; Patel, K.; Brosda, S.; et al. Identification of the CIMP-like subtype and aberrant methylation of members of the chromosomal segregation and spindle assembly pathways in esophageal adenocarcinoma. Carcinogenesis 2016, 37, 356–365. [Google Scholar] [CrossRef] [Green Version]
  38. Zhou, W.; Triche, T.J.; Laird, P.W.; Shen, H. SeSAMe: Reducing artifactual detection of DNA methylation by Infinium BeadChips in genomic deletions. Nucleic Acids Res. 2018, 46, e123. [Google Scholar] [CrossRef] [Green Version]
  39. Di Lena, P.; Sala, C.; Prodi, A.; Nardini, C. Missing value estimation methods for DNA methylation data. Bioinformatics 2019, 35, 3786–3793. [Google Scholar] [CrossRef]
  40. Morris, T.J.; Butcher, L.M.; Feber, A.; Teschendorff, A.E.; Chakravarthy, A.R.; Wojdacz, T.K.; Beck, S. ChAMP: 450k Chip Analysis Methylation Pipeline. Bioinformatics 2014, 30, 428–430. [Google Scholar] [CrossRef] [Green Version]
  41. Hoadley, K.A.; Yau, C.; Wolf, D.M.; Cherniack, A.D.; Tamborero, D.; Ng, S.; Leiserson, M.D.M.; Niu, B.; McLellan, M.D.; Uzunangelov, V.; et al. Multiplatform Analysis of 12 Cancer Types Reveals Molecular Classification within and across Tissues of Origin. Cell 2014, 158, 929–944. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Sanchez-Vega, F.; Mina, M.; Armenia, J.; Chatila, W.K.; Luna, A.; La, K.C.; Dimitriadoy, S.; Liu, D.L.; Kantheti, H.S.; Saghafinia, S.; et al. Oncogenic Signaling Pathways in The Cancer Genome Atlas. Cell 2018, 173, 321–337.e310. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. He, Y.; Zhou, H.; Xu, H.; You, H.; Cheng, H. Construction of an Immune-Related lncRNA Signature That Predicts Prognosis and Immune Microenvironment in Osteosarcoma Patients. Front. Oncol. 2022, 12, 769202. [Google Scholar] [CrossRef] [PubMed]
  44. Xue, W.; Chen, J.; Liu, X.; Gong, W.; Zheng, J.; Guo, X.; Liu, Y.; Liu, L.; Ma, J.; Wang, P.; et al. PVT1 regulates the malignant behaviors of human glioma cells by targeting miR-190a-5p and miR-488-3p. Biochim. Biophys. Acta (BBA) Mol. Basis Dis. 2018, 1864, 1783–1794. [Google Scholar] [CrossRef]
  45. Xu, S.; Wang, T.; Song, W.; Jiang, T.; Zhang, F.; Yin, Y.; Jiang, S.-W.; Wu, K.; Yu, Z.; Wang, C.; et al. The inhibitory effects of AR/miR-190a/YB-1 negative feedback loop on prostate cancer and underlying mechanism. Sci. Rep. 2015, 5, 13528. [Google Scholar] [CrossRef] [Green Version]
  46. Yu, Y.; Cao, X.-C. miR-190-5p in human diseases. Cancer Cell Int. 2019, 19, 257. [Google Scholar] [CrossRef]
  47. Li, H.; Zhang, X.; Jin, Z.; Yin, T.; Duan, C.; Sun, J.; Xiong, R.; Li, Z. MiR-122 Promotes the Development of Colon Cancer by Targeting ALDOA In Vitro. Technol. Cancer Res. Treat. 2019, 18, 1533033819871300. [Google Scholar] [CrossRef] [Green Version]
  48. Lin, X.; Han, L.; Gu, C.; Lai, Y.; Lai, Q.; Li, Q.; He, C.; Meng, Y.; Pan, L.; Liu, S.; et al. MiR-452-5p promotes colorectal cancer progression by regulating an ERK/MAPK positive feedback loop. Aging 2021, 13, 7608–7626. [Google Scholar] [CrossRef]
  49. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  50. Singla, A.; Wang, J.; Yang, R.; Geller, D.S.; Loeb, D.M.; Hoang, B.H. Wnt Signaling in Osteosarcoma. Adv. Exp. Med. Biol. 2020, 8, 125–139. [Google Scholar] [CrossRef]
  51. Hensel, J.; Wetterwald, A.; Temanni, R.; Keller, I.; Riether, C.; van der Pluijm, G.; Cecchini, M.G.; Thalmann, G.N. Osteolytic cancer cells induce vascular/axon guidance processes in the bone/bone marrow stroma. Oncotarget 2018, 9, 28877–28896. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Ben-Ghedalia-Peled, N.; Vago, R. Wnt Signaling in the Development of Bone Metastasis. Cells 2022, 11, 3934. [Google Scholar] [CrossRef]
  53. Jalanko, T.; de Jong, J.J.; Gibb, E.A.; Seiler, R.; Black, P.C. Genomic Subtyping in Bladder Cancer. Curr. Urol. Rep. 2020, 21, 9. [Google Scholar] [CrossRef] [PubMed]
  54. Collisson, E.A.; Bailey, P.; Chang, D.K.; Biankin, A.V. Molecular subtypes of pancreatic cancer. Nat. Rev. Gastroenterol. Hepatol. 2019, 16, 207–220. [Google Scholar] [CrossRef]
  55. Song, Y.-J.; Xu, Y.; Deng, C.; Zhu, X.; Fu, J.; Chen, H.; Lu, J.; Xu, H.; Song, G.; Tang, Q.; et al. Gene Expression Classifier Reveals Prognostic Osteosarcoma Microenvironment Molecular Subtypes. Front. Immunol. 2021, 12, 623762. [Google Scholar] [CrossRef] [PubMed]
  56. Zhang, Y.; Chen, F.; Creighton, C.J. Pan-cancer molecular subtypes of metastasis reveal distinct and evolving transcriptional programs. Cell Rep. Med. 2023, 4, 100932. [Google Scholar] [CrossRef]
  57. Kong, C.; Hansen, M.F. Biomarkers in osteosarcoma. Expert Opin. Med Diagn. 2009, 3, 13–23. [Google Scholar] [CrossRef] [Green Version]
  58. Tirabosco, R.; Berisha, F.; Ye, H.; Halai, D.; Amary, M.F.; Flanagan, A.M. Assessment of MUC4 expression in primary bone tumours. Histopathology 2013, 63, 142–143. [Google Scholar] [CrossRef]
  59. Liu, W.; Huang, Z.; Ding, Y.; Yang, Y.; Jin, T.; Deng, Z.; Xu, H.; Wang, X.; Zhang, Y.; Wang, H.; et al. Comprehensive genomic profiling of patients with favorable prognosis in osteosarcoma. J. Clin. Oncol. 2019, 37, e22501. [Google Scholar] [CrossRef]
  60. Chen, T.; Chen, Z.; Lian, X.; Wu, W.; Chu, L.; Zhang, S.; Wang, L. MUC 15 Promotes Osteosarcoma Cell Proliferation, Migration and Invasion through Livin, MMP-2/MMP-9 and Wnt/β-Catenin Signal Pathway. J. Cancer 2021, 12, 467–473. [Google Scholar] [CrossRef]
  61. Chiappetta, C.; Puggioni, C.; Carletti, R.; Petrozza, V.; Della Rocca, C.; Di Cristofano, C. The nuclear-cytoplasmic trafficking of a chromatin-modifying and remodelling protein (KMT2C), in osteosarcoma. Oncotarget 2018, 9, 30624–30634. [Google Scholar] [CrossRef] [Green Version]
  62. Chiappetta, C.; Mancini, M.; Lessi, F.; Aretini, P.; De Gregorio, V.; Puggioni, C.; Carletti, R.; Petrozza, V.; Civita, P.; Franceschi, S.; et al. Whole-exome analysis in osteosarcoma to identify a personalized therapy. Oncotarget 2017, 8, 80416–80428. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Du, X.; Yang, J.; Yang, D.; Tian, W.; Zhu, Z. The genetic basis for inactivation of Wnt pathway in human osteosarcoma. BMC Cancer 2014, 14, 450. [Google Scholar] [CrossRef]
  64. Matsuoka, K.; Bakiri, L.; Wolff, L.I.; Linder, M.; Mikels-Vigdal, A.; Patiño-García, A.; Lecanda, F.; Hartmann, C.; Sibilia, M.; Wagner, E.F. Wnt signaling and Loxl2 promote aggressive osteosarcoma. Cell Res. 2020, 30, 885–901. [Google Scholar] [CrossRef]
  65. Wang, L.; Jin, F.; Qin, A.; Hao, Y.; Dong, Y.; Ge, S.; Dai, K. Targeting Notch1 signaling pathway positively affects the sensitivity of osteosarcoma to cisplatin by regulating the expression and/or activity of Caspase family. Mol. Cancer 2014, 13, 139. [Google Scholar] [CrossRef] [Green Version]
  66. McManus, M.M.; Weiss, K.R.; Hughes, D.P.M. Understanding the Role of Notch in Osteosarcoma. Adv. Exp. Med. Biol. 2014, 804, 67–92. [Google Scholar] [CrossRef]
  67. Wang, D.-Y.; Wu, Y.-N.; Huang, J.-Q.; Wang, W.; Xu, M.; Jia, J.-P.; Han, G.; Mao, B.-B.; Bi, W.-Z. Hippo/YAP signaling pathway is involved in osteosarcoma chemoresistance. Chin. J. Cancer 2016, 35, 47. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  68. Fang, F.; VanCleave, A.; Helmuth, R.; Torres, H.; Rickel, K.; Wollenzien, H.; Sun, H.; Zeng, E.; Zhao, J.; Tao, J. Targeting the Wnt/β-catenin pathway in human osteosarcoma cells. Oncotarget 2018, 9, 36780–36792. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  69. Dallmayer, M.; Li, J.; Ohmura, S.; Rubio, R.A.; Baldauf, M.C.; Hölting, T.L.B.; Musa, J.; Knott, M.M.L.; Stein, S.; Cidre-Aranaz, F.; et al. Targeting the CALCB/RAMP1 axis inhibits growth of Ewing sarcoma. Cell Death Dis. 2019, 10, 116. [Google Scholar] [CrossRef] [Green Version]
  70. Logan, M.; Anderson, P.D.; Saab, S.T.; Hameed, O.; Abdulkadir, S.A. RAMP1 Is a Direct NKX3.1 Target Gene Up-Regulated in Prostate Cancer that Promotes Tumorigenesis. Am. J. Pathol. 2013, 183, 951–963. [Google Scholar] [CrossRef] [Green Version]
  71. Fu, P.; Shi, Y.; Chen, G.; Gao, Z. Establishment of key genes and associated outcomes in osteosarcoma patients using bioinformatics methods. Int. J. Clin. Exp. Med. 2021, 14, 2242–2249. [Google Scholar]
  72. Shi, D.; Mu, S.; Pu, F.; Liu, J.; Zhong, B.; Hu, B.; Ni, N.; Wang, H.; Luu, H.H.; Haydon, R.C.; et al. Integrative analysis of immune-related multi-omics profiles identifies distinct prognosis and tumor microenvironment patterns in osteosarcoma. Mol. Oncol. 2022, 16, 2174–2194. [Google Scholar] [CrossRef] [PubMed]
  73. Baumhoer, D.; Elsner, M.; Smida, J.; Zillmer, S.; Rauser, S.; Schoene, C.; Balluff, B.; Bielack, S.; Jundt, G.; Walch, A.; et al. CRIP1 expression is correlated with a favorable outcome and less metastases in osteosarcoma patients. Oncotarget 2011, 2, 970–975. [Google Scholar] [CrossRef] [Green Version]
  74. Cassoni, P.; Allia, E.; Marrocco, T.; Ghè, C.; Ghigo, E.; Muccioli, G.; Papotti, M. Ghrelin and cortistatin in lung cancer: Expression of peptides and related receptors in human primary tumors and in vitro effect on the H345 small cell carcinoma cell line. J. Endocrinol. Investig. 2006, 29, 781–790. [Google Scholar] [CrossRef] [PubMed]
  75. Cassoni, P.; Muccioli, G.; Marrocco, T.; Volante, M.; Allia, E.; Ghigo, E.; Deghenghi, R.; Papotti, M. Cortistatin-14 inhibits cell proliferation of human thyroid carcinoma cell lines of both follicular and parafollicular origin. J. Endocrinol. Investig. 2002, 25, 362–368. [Google Scholar] [CrossRef]
  76. Wu, Z.-L.; Deng, Y.-J.; Zhang, G.-Z.; Ren, E.-H.; Yuan, W.-H.; Xie, Q.-Q. Development of a novel immune-related genes prognostic signature for osteosarcoma. Sci. Rep. 2020, 10, 18402. [Google Scholar] [CrossRef]
Cancers 15 02134 g001 550
Figure 1. (A) Oncoplot showing the top 5 altered genes with more than 10% samples containing mutations in the TCGA enriched pathways. The fraction and the number of samples with variations in the corresponding genes are shown on the right as percentages and the bar plots (B) Enriched oncogenic pathways showing the fraction of genes and samples affected.
Figure 1. (A) Oncoplot showing the top 5 altered genes with more than 10% samples containing mutations in the TCGA enriched pathways. The fraction and the number of samples with variations in the corresponding genes are shown on the right as percentages and the bar plots (B) Enriched oncogenic pathways showing the fraction of genes and samples affected.
Cancers 15 02134 g001
Cancers 15 02134 g002 550
Figure 2. Kaplan-Meier plot for high vs low expression group in TARGET OS data with p-value from log-rank test and Cox regression model for genes, METTL20, MYH10, and PDE1B with HR < 1 (A), and genes, RAMP1, TAC4, and MT1 with HR > 1 (B).
Figure 2. Kaplan-Meier plot for high vs low expression group in TARGET OS data with p-value from log-rank test and Cox regression model for genes, METTL20, MYH10, and PDE1B with HR < 1 (A), and genes, RAMP1, TAC4, and MT1 with HR > 1 (B).
Cancers 15 02134 g002
Cancers 15 02134 g003 550
Figure 3. (A) Correlation plot showing significantly correlated canonical pathway genes (p < 0.05). Negative correlated genes are shown in red, and positive in blue. (B) List of 34 genes found significantly associated with survival (p < 0.05), their associated canonical pathways, and corresponding HR values.
Figure 3. (A) Correlation plot showing significantly correlated canonical pathway genes (p < 0.05). Negative correlated genes are shown in red, and positive in blue. (B) List of 34 genes found significantly associated with survival (p < 0.05), their associated canonical pathways, and corresponding HR values.
Cancers 15 02134 g003
Cancers 15 02134 g004 550
Figure 4. Average silhouette width and NMF classification with rank using 190 miRNAs, 1107 mutations, 2500 protein-coding regions, 500 lncRNA, and 5000 CpG methylation probe data for 85 OS-TARGET patients which was further used as input for COCA for obtaining the final two molecular subtypes.
Figure 4. Average silhouette width and NMF classification with rank using 190 miRNAs, 1107 mutations, 2500 protein-coding regions, 500 lncRNA, and 5000 CpG methylation probe data for 85 OS-TARGET patients which was further used as input for COCA for obtaining the final two molecular subtypes.
Cancers 15 02134 g004
Cancers 15 02134 g005 550
Figure 5. (A) Volcano plot showing the differentially expressed genes between the two molecular subtypes. The downregulated (green) and upregulated (red) genes obtained at abs(log2Fold Change) > 1.5 and adjusted p-value < 0.05 are highlighted. (B) Heatmap showing top 100 most variable differentially expressed genes obtained at abs(log2Fold Change) > 1 and adjusted p-value < 0.01 cutoff found between the two molecular clusters.
Figure 5. (A) Volcano plot showing the differentially expressed genes between the two molecular subtypes. The downregulated (green) and upregulated (red) genes obtained at abs(log2Fold Change) > 1.5 and adjusted p-value < 0.05 are highlighted. (B) Heatmap showing top 100 most variable differentially expressed genes obtained at abs(log2Fold Change) > 1 and adjusted p-value < 0.01 cutoff found between the two molecular clusters.
Cancers 15 02134 g005
Table
Table 1. List of the top 49 significant (p-value < 0.001) genes found in survival analysis. The gene list is sorted based on beta values. HR: Hazard Ratio.
Table 1. List of the top 49 significant (p-value < 0.001) genes found in survival analysis. The gene list is sorted based on beta values. HR: Hazard Ratio.
GeneHRHR_low95HR_up95BetaGeneHRHR_low95HR_up95Beta
RAMP15.72.3141.7LY860.240.10.56−1.4
TAC45.52.6121.7SHISA50.240.10.57−1.4
TIMM23B4.72.1111.6FOLR20.250.110.56−1.4
COL5A24.72111.6C11orf450.250.110.57−1.4
MT1A4.82.1111.6UBE2L30.240.110.54−1.4
CRIP14.81.9121.6CYFIP10.240.10.55−1.4
UFC14.61.9111.5SNX10.220.0910.55−1.5
PROSER24.41.9101.5ACTB0.220.0950.5−1.5
NHEJ14.72111.5MPP10.220.0910.55−1.5
CORT4.62.1101.5CD1800.230.0920.56−1.5
CKMT24.128.81.4TPMT0.220.090.55−1.5
CHD1L3.91.78.91.4BBS40.220.0880.53−1.5
LGR64.11.99.11.4ITGAM0.220.0910.56−1.5
IFFO24.21.89.71.4TCN20.230.0930.56−1.5
MAFK4.11.89.71.4SIRPA0.210.0870.53−1.5
PGAM43.91.88.51.4SETD90.190.0790.48−1.6
SNAP913.81.881.3SNTB20.210.0850.53−1.6
PDE4C3.51.77.51.3MITF0.210.0790.55−1.6
KRT23.61.77.61.3NUBP10.170.0730.41−1.8
GCSAM3.71.781.3IRF2BPL0.170.0690.42−1.8
BAI13.71.77.81.3ERCC40.170.0660.46−1.8
CFAP443.81.881.3PDE1B0.150.0450.49−1.9
FAM166B3.51.77.31.3MYH100.160.0540.46−1.9
APBB1IP0.240.110.55−1.4METTL200.140.0530.36−2
COMMD90.240.110.55−1.4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Southekal, S.; Shakyawar, S.K.; Bajpai, P.; Elkholy, A.; Manne, U.; Mishra, N.K.; Guda, C. Molecular Subtyping and Survival Analysis of Osteosarcoma Reveals Prognostic Biomarkers and Key Canonical Pathways. Cancers 2023, 15, 2134. https://0-doi-org.brum.beds.ac.uk/10.3390/cancers15072134

AMA Style

Southekal S, Shakyawar SK, Bajpai P, Elkholy A, Manne U, Mishra NK, Guda C. Molecular Subtyping and Survival Analysis of Osteosarcoma Reveals Prognostic Biomarkers and Key Canonical Pathways. Cancers. 2023; 15(7):2134. https://0-doi-org.brum.beds.ac.uk/10.3390/cancers15072134

Chicago/Turabian Style

Southekal, Siddesh, Sushil Kumar Shakyawar, Prachi Bajpai, Amr Elkholy, Upender Manne, Nitish Kumar Mishra, and Chittibabu Guda. 2023. "Molecular Subtyping and Survival Analysis of Osteosarcoma Reveals Prognostic Biomarkers and Key Canonical Pathways" Cancers 15, no. 7: 2134. https://0-doi-org.brum.beds.ac.uk/10.3390/cancers15072134

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop