Deep RNA sequencing reveals the dynamic regulation of miRNA, lncRNAs, and mRNAs in osteosarcoma tumorigenesis and pulmonary metastasis

Osteosarcoma (OS) is the most common pediatric malignant bone tumor, and occurrence of pulmonary metastasis generally causes a rapid and fatal outcome. Here we aimed to provide clues for exploring the mechanism of tumorigenesis and pulmonary metastasis for OS by comprehensive analysis of microRNA (miRNA), long non-coding RNA (lncRNA), and mRNA expression in primary OS and OS pulmonary metastasis. In this study, deep sequencing with samples from primary OS (n = 3), pulmonary metastatic OS (n = 3), and normal controls (n = 3) was conducted and differentially expressed miRNAs (DEmiRNAs), lncRNAs (DElncRNAs), and mRNAs (DEmRNAs) between primary OS and normal controls as well as pulmonary metastatic and primary OS were identified. A total of 65 DEmiRNAs, 233 DElncRNAs, and 1405 DEmRNAs were obtained between primary OS and normal controls; 48 DEmiRNAs, 50 DElncRNAs, and 307 DEmRNAs were obtained between pulmonary metastatic and primary OS. Then, the target DEmRNAs and DElncRNAs regulated by the same DEmiRNAs were searched and the OS tumorigenesis-related and OS pulmonary metastasis-related competing endogenous RNA (ceRNA) networks were constructed, respectively. Based on these ceRNA networks and Venn diagram analysis, we obtained 3 DEmiRNAs, 15 DElncRNAs, and 100 DEmRNAs, and eight target pairs including miR-223-5p/(CLSTN2, AC009951.1, LINC01705, AC090673.1), miR-378b/(ALX4, IGSF3, SULF1), and miR-323b-3p/TGFBR3 were involved in both tumorigenesis and pulmonary metastasis of OS. The TGF-β superfamily co-receptor TGFBR3, which is regulated by miR-323b-3p, acts as a tumor suppressor in OS tumorigenesis and acts as a tumor promoter in pulmonary metastatic OS via activation of the epithelial–mesenchymal transition (EMT) program. In conclusion, the OS transcriptome (miRNA, lncRNA, and mRNA) is dynamically regulated. These analyses might provide new clues to uncover the molecular mechanisms and signaling networks that contribute to OS progression, toward patient-tailored and novel-targeted treatments.


Introduction
Osteosarcoma (OS) is one of the main primary malignant bone tumor subtypes which mostly occurs in adolescents at sites of rapid bone growth 1 . Although intensive efforts to improve both chemotherapeutics and surgical management have been made, the high local aggressiveness and rapid metastasizing potential to lung results in poor survival for patients with OS. Therefore, the ultimate treatment depends on primary OS control and the removal of small metastases. OS is a pathology that affects bone remodeling, involving alterations in both osteoblast and osteoclast functions. However, the mechanisms underlying its initiation and progression remain unclear.
Non-coding RNAs (ncRNAs) have no ability of coding proteins while they can act as functional RNAs. Based on the transcript size, ncRNAs are grouped into small ncRNAs (<200 bp) and long ncRNAs (>200 bp, up to 100 kb). MicroRNAs (miRNAs), a class of small ncRNAs (≈22 nt), are crucial to the regulation of gene expression through partial base-pairing with target mRNAs. They have multiple roles in various biological processes that affect basic cellular functions, including cell proliferation, differentiation, death, and tumorigenesis 2 . Unlike miRNAs, long non-coding RNA (lncRNAs) play critical and complicated roles in the regulation of various biological processes, including chromatin modification, transcription, and post-transcriptional processing 3 . Currently, growing evidences indicated that there are interactions between lncRNAs and miRNAs, the downstream target genes of which have been closely related to tumor pathogenesis.
RNA sequencing (RNA-seq) has been used widely to study specific gene expression patterns at different developmental stages. In this study, we obtained the miRNA, mRNA, and lncRNA expression data from normal controls, primary OS, and pulmonary metastatic OS based on RNA-seq, and we constructed the competing endogenous (ceRNA) network to elaborate the interactions and potential crosstalk between the differentially expressed hub lncRNAs (DElncRNAs), miRNAs (DEmiRNAs), and mRNAs (DEmRNAs).

Sample preparation
In the present study, we recruited patients with OS from the Third Affiliated Hospital of Kunming Medical University. Detailed information of patients is displayed in Table 1. The fresh tumor tissues were obtained from the primary lesion of three patients with primary OS (P1-P3) and three patients with pulmonary metastatic OS (M1-M3) after surgical resection. The control noncancerous tissues were obtained from distal tumor location of three patients with primary OS (N1-N3). Tissue samples were frozen in liquid nitrogen and stored at −80°C before RNA isolation. The present study complied with Declaration of Helsinki and was approved by the Institutional Review Boards of the Third Affiliated Hospital of Kunming Medical University. Moreover, all subjects provided written informed consent.

RNA extraction and quality monitoring
All the surgical specimens were subjected to RNA extraction using the Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The RNA quality was evaluated with the NanoDrop2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies). Purified RNA was stored at −80°C until required.
Small RNA library construction, sequencing, and data processing Following extraction and purification, about 1 μg total RNA per sample was used to construct the small RNA (sRNA) library using TrueSeq small RNA library prep kit (Illumina San Diego CA, USA) according to the manufacturer's instruction. Adapters were ligated to the 3′ end of the RNA, followed by the ligation of the 5′ adapter. Subsequently, the RNA was reverse transcribed to create single-stranded cDNA, followed by single-end sequencing (50 base pairs in length) on an Illumina on the HiSeq4000 sequencer (Illumina, San Diego, CA, USA).
lncRNA + mRNA sequencing and data processing A total of 3 μg RNA per sample was used for the RNA sample preparations. After removing the ribosomal RNA, the rRNA-depleted RNA was fragmented and the cDNA library was constructed using the Truseq RNA sample Prep Kit (Illumina, Inc., San Diego, CA, USA). The libraries were sequenced on an Illumina Hiseq 2500 platform (Illumina Inc., San Diego, CA, USA) according to the manufacturer's instructions and 125 bp paired-end reads were generated.
Raw reads of fastq format were then processed through in-house perl scripts. After triming the raw reads, we obtained the clean reads, which were mapped to the human reference genome Ensembl V84 using Tophat. The mapped reads were quantified with cuffquant, and the differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) between samples were identified using Cuffdiff program from the Cufflinks package. The p-value < 0.01 and |log2 (Fold_change)| > 2 were used as the cut-off criteria.

Function enrichment analysis
We performed Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses on these DEmRNAs and predicted target genes of DEmiRNAs and DElncRNAs. GO term and KEGG pathway analyses of coding genes were performed using GeneCodis3 bioinformatics resources 6 . Both GO terms and KEGG pathways with corrected. p-Values <0.05 were considered to be significantly enriched.
Schema for integrative analysis of DEmiRNAs, DEncRNA, and DEmRNA Systematic bioinformatic analysis was developed based on possible functional relationships between DEmiRNAs, DEncRNA, and DEmRNA. Firstly, by scanning for conserved miRNA target sites with RNA22, miRanda, miRDB, miRWalk, PICTAR2, and Targetscan, we predicted the target genes and target lncRNAs for the DEmiRNAs. Secondly, we searched coding genes within the 100-kb upstream and downstream regions of each DElncRNAs and found the cis-acting genes. According to the functional relationships between these molecules, the miRNAtarget gene regulatory network, miRNA-lncRNA target regulatory network, lncRNA-mRNA co-expression network were established, respectively. Next, we constructed the ceRNA network.

Quantitative real-time polymerase chain reaction (qRT-PCR)
To validate the expression levels of the selected lncRNAs by qRT-PCR, RNA samples from the additional 36 individuals with primary OS and 33 individuals with pulmonary metastatic OS were collected. To validate the expression levels of selected genes and miRNAs by qRT-PCR, RNA samples from the additional 30 individuals with primary OS and 27 individuals with pulmonary metastatic OS were collected. Total RNA was isolated by using the Trizol reagent (Invitrogen, USA) according to manufacturer's protocol. The mRNA template was reversely transcribed into cDNA using reverse transcriptase Kit (TaKaRa, Dalian, China). The miRNA reverse transcription was performed using miRcute miRNA First-strand cDNA Synthesis kits (TIANGEN, China). Forward and reverse primers were designed and qRT-PCR was carried out on BIO-RAD IQ5 RT-PCR Detection System (Bio-Rad Laboratories Inc., Germany). The expression levels of the selected lncRNAs and miR-NAs were normalized against the snU6. The expression levels of the selected genes were normalized against GAPDH.

Cross-validation
The miRNA expression data of GSE65071 was downloaded from GEO database (https://www.ncbi.nlm.nih. gov/geo/), including 20 plasma samples from OS cases and 15 plasma samples from controls plasma. The DEmiRNAs were validated between comparison of case group and control samples. Moreover, the mRNA data of GSE14359 was also downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/), including 10 conventional OS tissues and eight OS lung metastasis tissues. The DEmRNAs were validated between comparison of OS lung metastasis group and OS group.

Sequencing and mapping of the OS transcriptome
We sequenced the cDNA and sRNA libraries of nine tissue samples from three primary OS patients, three pulmonary metastatic OS patients, and three controls. Counts of clean reads and mapped ratio of sequencing results were displayed in Supplemental Table S1. The overall workflow is shown in Fig. 1.

Deep RNA-seq revealed distinct expression signatures of coding and ncRNAs in OS progression
Principal component analyses (PCA) revealed that miRNA, lncRNA, and mRNA expression profiles distinguish primary OS, pulmonary metastatic OS from the controls (Fig. 2a-c). Using the criterion of p < 0.01 and |log2(fold change)| > 2, we detected 65 DEmiRNAs, 233 DElncRNAs, and 1405 DEmRNAs in primary OS compared with the normal controls. Top ten miRNAs, lncRNAs, and mRNAs exhibiting significant up-and downregulation are listed in Table 2. Totally, we detected 48 DEmiRNAs, 50 DElncRNAs, and 307 DEmRNAs in pulmonary metastatic OS compared with primary OS. Top ten miRNAs, lncRNAs, and mRNAs exhibiting significant up-and downregulation are listed in Table 2. Unsupervised hierarchical clustering of the DEmiRNAs, DElncRNAs, and DEmRNAs ( Fig. 2d-f) revealed a distinct expression signature of all three RNA species in primary OS and OS pulmonary metastasis, compared to the control samples.
The distribution of DElncRNAs and DEmRNAs in primary OS compared with the normal controls is illustrated in Fig. 2g. The unsupervised clustering showed two robust clusters: one cluster encompassing all of primary OS and another cluster containing all of the controls. This indicated that tumors and controls might have different expression patterns (Fig. 2h-j). Moreover, the chromosomes distribution of DElncRNAs and DEmRNAs in pulmonary metastatic OS compared with primary OS is illustrated in Fig. 2k. Unsupervised hierarchical clustering of the expression profiles of mRNAs, miRNAs, and lncRNAs ( Fig. 2l-n) revealed that lncRNA and mRNA expression profiles can largely distinguish the pulmonary metastatic OS and primary OS.
Previous studies have reported that lncRNAs may act in cis and affect the gene expression of their chromosomal neighborhood 45 , and most lncRNA transcripts can also be derived from divergent transcription 46 . We annotated the location relationship between each lncRNA and its cis target genes and obtained 125 lncRNAs and their neighboring genes pairs in total. Accordingly, the lncRNA-mRNA networks were constructed and visualized (Supplemental Fig. 3).
According to the target pairs of miRNA-mRNA, miRNA-lncRNA, and lncRNA-cis target gene, we constructed a ceRNA network (Fig. 3a). In particular, three miRNAs (miR-223-5p, miR-378b, and miR-323b-3p) were not only DEmiRNAs between primary OS and normal control but also DEmiRNAs between pulmonary metastasis OS and primary OS, which suggested that these three miRNAs might involve in both oncogenesis and metastasis of OS. Hence, subnetworks of miR-223-5p, miR-378b, and miR-323b-3p were shown in Fig. 3b. Functional annotation of mRNAs in the ceRNA network found that these mRNAs were mainly involved in seven GO terms (Fig. 3c), and KEGG pathways including Phagosome, ECM-receptor interaction, and two well-established cancer pathways, apoptosis 47 and MAPK signaling pathway 48 (Fig. 3d and Supplementary Table S2).
Totally, 43 lncRNAs, 43 miRNAs, and 143 mRNAs were involved in the ceRNA network (Fig. 4a). Since Fig. 3 The OS tumorigenesis-related ceRNA network and their characteristics. a Global view of the ceRNA network. This network consists of 96 DElncRNAs, 50 DEmiRNAs, and 125 DEmRNAs. b Subnetworks of miR-223-5p, miR-378b, and miR-323b-3p. In Fig. 4a, b, rectangles, diamonds, and ellipses represented DEmiRNAs, DElncRNAs, and DEmRNAs between primary OS and normal controls, respectively. Red and green color represented upregulation and downregulation in primary OS compared to normal controls. c The functional enrichment map of GO terms. d Significantly enriched KEGG pathway of mRNAs in the ceRNA network miR-223-5p, miR-378b, and miR-323b-3p were differentially expressed between primary OS and normal controls as well as pulmonary metastasis OS and primary OS. Subnetworks of miR-223-5p, miR-378b, and miR-323b-3p are shown in Fig. 4b. Functional annotation revealed that mRNAs in the ceRNA network were mainly involved in seven GO terms (Fig. 4c), and KEGG pathways including Hematopoietic cell lineage, Gap junction, Focal adhesion, Complement and coagulation cascades and a well-established cancer pathway, p53 signaling pathway 64 (Fig. 4d and Supplementary  Table S2).

qRT-PCR and cross-validation
Twelve lncRNAs randomly selected from top 15 DElncRNAs between pulmonary metastatic and primary OS were validated by qRT-PCR. Relative expression of DElncRNAs was determined by quantitative RT-PCR of 36 primary OS samples and 33 pulmonary metastatic OS samples (Fig. 6). Expression of these 12 DElncRNAs was consistent with that in RNA-seq results in this study. miRNA was validated using GSE65071 from GEO database (20 plasma from OS cases and 15 controls plasma). DEmiRNA (miR-223-5p) validated by GSE65071 were randomly selected from three shared DEmiRNAs (miR-223-5p, miR-378b, and miR-323b-3p) in pulmonary metastatic OS vs. primary OS and primary OS vs. normal controls (Fig. 6). Expression of miR-223-5p was consistent with that in RNA-seq results in this study.
Genes were validated by GSE14359 from GEO database (10 conventional OS tissues and eight OS lung metastasis tissues). Five genes validated by GSE14359 were randomly selected from targets of miR-223-5p, miR-378b, and miR-323b-3p which were DEmRNAs in pulmonary metastatic OS vs. primary OS and primary OS vs. normal controls (Fig. 6). Expression of these five genes was consistent with that in RNA-seq results in this study.
The central role of TGF-β superfamily co-receptor TGFBR3 in OS tumorigenesis and pulmonary metastatic OS Transforming growth factor-β (TGF-β) plays critical roles in the vicious cycle between OS cells and the bone tumor microenvironment, thus contributing to tumor development and lung metastases dissemination 65 . Our RNA-seq results showed that TGFBR3, a TGF-β superfamily co-receptor, was decreased in the process of OS tumorigenesis and increased in the process of OS pulmonary metastasis. In addition, TGFBR3 was a target of miR-323b-3p, and miR-323b-3p was increased in OS development while decreased in pulmonary metastatic OS based on our RNA-seq results. Moreover, both upregulated TGFBR3 and downregulated miR-323b-3p were observed in pulmonary metastatic OS compared to primary OS in our qRT-PCR results (Fig. 7). We established the central role of TGFBR3 in the OS development, which was shown in Fig. 7. At the early stage, miR-323b-3p inhibits the expression of TGFBR3, which acts as tumor suppressors in OS tumorigenesis by decreasing the TGF-β signaling. However, at the malignant stage, miR-323b-3p promotes the expression of TGFBR3, which acts as tumor promoters in pulmonary metastatic OS by enhancing the TGF-β signaling (Fig. 7). Previous study revealed that TGF-βs may act as a tumor suppressor by inhibiting the proliferation of epithelial cells and act as tumor promoters during the late stages of carcinogenesis by inducing epithelial-mesenchymal transition (EMT), to stimulate angiogenesis, and to favor immune evasion 65 . Hence, expression of EMT-promoting transcription factors, including ZEB1, ZEB2, TWIST1, and SNAI2 66 , was validated by qRT-PCR. All these four EMT-promoting transcription factors were elevated in pulmonary metastatic OS compared to primary OS based on qRT-PCR results (Fig. 7).

Discussion
OS is the most common pediatric malignant bone tumor with early pulmonary metastasis formation as a frequent occurrence. Once metastasized to the lung, OS generally causes a rapid and fatal outcome. To improve this situation, increasing attention has been given to identify the exact regulatory mechanism of OS development and malignancy. Recent years, ncRNAs have been found to be associated with wide range of biological regulatory functions 67 . The present study utilized nextgeneration sequencing to provide a quantitative and comprehensive analysis of the coding and non-coding transcriptome in primary OS and pulmonary metastasis. These analyses revealed significant differences in the patterns of miRNA, lncRNA, and mRNA expression in primary OS and pulmonary metastasis, as well as the dynamic changes of DEmiRNA, DElncRNA, and DEmRNA. Here we showed for the first time that the expression patterns of lncRNAs and mRNAs are more suitable to discriminate the controls, primary OS, and pulmonary metastatic OS samples.
In general, our data suggest that distinct populations of miRNAs, lncRNAs, and mRNAs are involved in the pathogenesis of primary OS and OS pulmonary Definite expression of key miRNAs, lncRNAs, and mRNAs and their dynamic expression trends. The definite expression was shown based on log 10 (FPKM) or log 10 (normalized read counts) value. g-i Circos plots of the pathways predicted to be targeted by miR-223-5p (g), miR-323b-3p (h), and miR-378b (i) metastasis. Based on the RNA sequence data, we identified 65 miRNAs, 233 lncRNAs, and 1405 mRNAs were differentially expressed in primary OS compared with the normal controls. We found that most of the DEmiRNAs may be associated with OS tumorigenesis. However, most of the DElncRNAs were with unknown function, which is mainly due to the few researches for them. Based on tumorigenesis-related ceRNAs, Apoptosis and MAPK Fig. 7 The central role of TGFBR3 in the OS development. TGFBR3 acts as a tumor suppressor in OS tumorigenesis, whereas TGFBR3 acts as a tumor promoter in OS pulmonary metastasis. Expressions of TGFBR3, miR-323b-3p, and four EMT-promoting transcription factors (ZEB1, ZEB2, TWIST1, and SNAI2) have been validated by qRT-PCR. Relative expression was determined by quantitative RT-PCR of 30 primary OS samples and 27 pulmonary metastatic OS samples signaling pathway were two significantly enriched pathways which were well-established cancer-related pathways 47,48 . PPP3CC, protein phosphatase 3 catalytic subunit gamma was a shared gene in these two pathways, decreased PPP3CC has been found in prostate cancer and gliomas 68,69 . We firstly found downregulation of PPP3CC in primary OS compared to normal controls in this study. Moreover, PPP3CC was a neighboring gene of a downregulated lncRNA in primary OS compared to normal controls, AC037459.2 (ENSG00000251034). PPP3CC-AC037459.2 interaction was speculated to involve with the processes of OS and other cancers by regulating Apoptosis and MAPK signaling pathway. The precise role of PPP3CC-AC037459.2 interaction in cancers needs further research. Besides, dysregulated mRNAs in the OS tumorigenesis-related ceRNA network were mainly involved in another two pathways, Phagosome and ECMreceptor interaction that highlighted their importance in OS. Recent studies identified that miRNA interactions with lncRNA and mRNA might play important roles in OS formation, pulmonary metastasis and prognosis, such as miR-30a 70 , miR-136 71 , miR-206 10 , miR-181b 32 , miR-29b-3p 72 , miR-29a 73 , miR-133a 74 , miR-224 30 , and miR-223 24 . In our ceRNA network, we also found these key miRNAs have most target lncRNAs or mRNAs. Therefore, our results suggested that these key miRNAs may play an important role in the progression and development of OS and the cancer genes related pathways.
Subsequently, we detected 48 miRNAs, 50 lncRNAs, and 307 mRNAs that had different expression patterns in pulmonary metastatic OS compared with primary OS. DEmRNAs in OS pulmonary metastasis-related ceRNA network were significantly enriched in a well-known cancer-related pathway, p53 signaling pathway 64 . SESN2, sestrin 2 was enriched in this pathway which was reported to involve with various cancers such as bladder, breast, and lung cancers [75][76][77] . In this study, SESN2 was firstly found to be upregulated in pulmonary metastatic OS compared to primary OS. These findings suggested that SESN2 might be a potential tumor suppressor of OS by regulating p53 signaling pathway.
Among which, the top ten miRNAs were considered as the most important ones participating in OS pulmonary metastasis. Previous studies also indicated that their dysregulation may contribute to the progression or OS metastasis, such as miR-20b 58 , miR-182 60 , miR-134 78 , and miR-183 79 . Our results first suggested that miR-612, miR-1197, miR-193b-3p, miR-1262, miR-144-3p, and miR-1269a may also play roles in OS metastasis.
Moreover, a total of three DEmiRNAs were inferred as the most promising candidate genes affecting OS development, which were further described as follows. Namløs et al. found that miR-223 was identified with an intermediate expression level in OS clinical samples compared to osteoblasts and bone 80 . MiR-223 may have a tumor suppressor function in OS through the PI3K/Akt/ mTOR pathway and could be used in anticancer therapies in OS 81 . miR-223/Ect2/p21 signaling is also an important pathway that regulates the OS cell cycle progression and proliferaion 82 . Combination of miR-223 downregulation and Ect2 upregulation may be a possible marker of poor prognosis in OS malignancy 24,83 . In our present study, we found that miR-223-5p was upregulated in primary OS, whereas downregulated in OS pulmonary metastasis. The results suggested that miR-223-5p may act as important roles in OS development. Grilli et al. observed the modulation of miR-378 using an OS differentiative model 15 . Novello et al. also found that miR-378 was significantly downregulated in OS vs. control, high-grade OS vs. lowgrade OS, and metastatic OS vs. non-metastatic OS patients by RT-PCR 84 . Our data showed that miR-378b was downregulated in primary OS, whereas upregulated in OS pulmonary metastasis. These observations also suggested that miR-378b may be essential in OS progression. Although the relationship of miR-323b-3p and OS was not reported, our study showed a significant dysregulation of miR-323b-3p during the OS development and malignancy, which implied its critical roles in OS.
TGF-β signaling pathway is critical in OS development and in their metastatic progression 65 . TGF-βs act as both tumor suppressors and tumor promoters, depending on the cancer type and tumor development timing 85 . Previous study revealed that TGF-βs may act as a tumor suppressor by inhibiting the proliferation of epithelial cells and act as tumor promoters during the late stages of carcinogenesis by inducing EMT, to stimulate angiogenesis, and to favor immune evasion 65 . Moreover, the increase of TGF-βs is also associated with the presence of metastases in lung 86 and is correlated with high-grade OS 87 . TGFBR3 is a co-receptor for the TGF-β superfamily, which can present ligand to the TGF-β signaling receptors. TGFBR3 is a tumor suppressor in many tissue types 88 . Our RNAseq results found that TGFBR3 is decreased in OS development and increased in pulmonary metastatic OS which was consistent with that in both the crossvalidation and qRT-PCR validation. Moreover, upregulated expression of several EMT-promoting transcription factors including ZEB1, ZEB2, TWIST1, and SNAI2 was found in pulmonary metastasis compared to primary OS based on our qRT-PCR validation results. We suggested that TGFBR3 acts as a tumor suppressor during the early stage of OS development and becomes a tumor promoter during the late stages of metastases via activation of the EMT program. Therefore, blocking TGF-β signaling may represent a promising therapeutic approach to treat OS patients. Additionally, TGFBR3 was a target of miR-323b-3p and miR-323b-3p was found to be upregulated in OS development while downregulated in pulmonary metastatic OS based on RNA-seq, cross-validation, and qRT-PCR validation results in this study. We speculated that TGF-β signaling pathway was regulated by miR-323b-3p in oncogenesis and metastasis of OS.
Our study has some limitations. The number of samples analyzed here was relatively small, and the samples were obtained from a heterogeneous cohort of patients and donors. It may introduce some bias. Moreover, little is known about the alteration and functional significance of lncRNAs and additional studies are needed to explore the functional roles of lncRNAs in OS. Taken together, our study revealed distinct relative abundance and expression pattern of miRNAs, lncRNAs, and mRNAs in human OS, highlighting the different biological roles of the individual RNA classes during OS progression. Our results provide valuable information for ncRNAs studies in the future.