Data collection
Lots of publicly available gene expression datasets have been submitted to GEO database (Gene Expression Omnibus Database). Here we carefully searched GEO database with the key words of “osteosarcoma AND mRNA AND microRNA”. Six datasets were screened out. Among the 6 datasets, three datasets (GSE89370, GSE89074 and GSE86109) only contain mRNA expression profiles no miRNA data, and two datasets (GSE89930 and GSE70415) contain both miRNA and mRNA expression profiles without detailed clinical characteristics. Only the dataset GSE39058 contain both miRNA and mRNA expression profiles and well documented clinical characteristics. This dataset was submitted by Katherine Hill in 2012 [34], which consists of 164 samples from 65 patients. The patients who showed < 80% tumor necrosis after chemotherapy were defined as chemoresistant, while who showed > = 80% tumor necrosis were defined as chemosensitive. Among the 164 samples, 37 samples contain mRNA expression profiles and 65 samples contain miRNA expression profiles. For mRNA expression profiles, 14 samples were chemosensitive and 23 were chemoresistant patients. For miRNA expression profiles, 32 sample were chemosensitive and 33 were chemoresistant. The experiment platforms for these data were Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip and Illumina Human v2 microRNA expression beadchip respectively.
Analysis of Differentially Expressed Genes
mRNA and miRNA expression profiles were divided into chemoresistant group and chemosensitive group. Then the expression profiles were processed using internal R scripts and available annotation databases. Briefly, background correction, quantile normalization and log2 transformation were firstly applied using GeneChip Robust Multi-array Analysis algorithm (GC-RMA) [35]. Then, useless probes were filtered out and mean expression values were calculated for each gene. Finally, differentially expressed genes (DEGs) and miRNAs were obtained by Limma package (Linear Models for Microarray Analysis) within Bioconductor [36]. Adjust p-value and absolute log2 fold change were set to less than or equal to 0.05 and greater or equal to 2 respectively.
Function Analysis
The differentially expressed genes and miRNAs along with fold change and p-value were uploaded to Ingenuity Pathway Analysis (IPA) (QIAGEN, Redwood City, CA, USA). Then core analysis was performed with parameters setting as direct and indirect relationships included, causal networks included, selected all node types, and include experimentally observed results. The rest parameters were set as default. The core analysis consists of canonical pathway enrichment, upstream regulator analysis, downstream functions and disease, regulator effect analysis, network analysis and so on. In the results, IPA will calculate a statistical z-score indicating effect direction, and usually z-score large than 0 means activation.
miRNA and mRNA Pairings
We have built mRNA-miRNA pairings by microRNA Target Filter module in IPA to construct the regulation network. The analysis is based on experimentally validated mRNA and miRNA interaction from TarBase [37], miRecords [38] and peer-reviewed literatures as well as predicted mRNA-miRNA interaction from TargetScan [39].
Survival analysis
To evaluate the prognostic value of identified genes, we downloaded the dataset of GSE21257 with 53 samples including documended clinical characteristics. The dataset was divided into two groups according to the quantile expression of the identified genes [40]. Survival analysis was carried out using the Kaplan–Meier method with the log-rank test to compare overall survival and gene expression groups.
Osteosarcoma tissue analysis
Comparing the pain degree, tumor size, tumor boundaries, and calcification before and after the chemotherapy was used to evaluate chemotherapy efficacy, which was classified into three categories: remission (RD), stability (SD), and progression (PD). RD was defined as pain relief, clear boundaries, no increase in tumor size and calcification. SD was defined as no increased pain, no change in tumor boundaries, and no increase in tumor size. PD was defined as increased pain, unclear tumor boundaries, increase in tumor size and calcification. RD represents chemotherapy sensitivity, while SD and PD represent drug resistance. Based on this criteria [41], the new enrolled 68 osteosarcoma patients were divided into chemosensitive (38 samples) and chemoresistant (30 samples) group. These samples were collected during May 2013 to Dec. 2018 from the department of Bone & Soft Tissue Surgery, Hubei Cancer Hospital. All patients enrolled provided consent and under approval by the Hubei Cancer Hospital Ethics Committee (ID: LLHBCH2018KY-018). All tissue samples were collected from patients during surgery or biopsy, cut into blocks with a diameter of 3–5 mm and transferred to liquid nitrogen within 30 min.
RT-PCR experiments were preformed on the 68 samples to validate the mRNA expression of the identified differentiately expressed genes (TIMP2 and BAX). Total RNA was extracted using TRIzol reagent with manufacturer’s instruction (ThermoFisher Inc., USA). Then cDNA was synthesized using M-MLV Reverse Transcriptase from Promega (Promega, Madison, WI, USA), and mRNA expression was detected using 7500 Real-Time PCR system (ThermoFisher, USA). The internal GAPDH mRNA expression was used for normalization, and relative quantification was calculated using 2−ΔCt method [42].
Cell lines and Lentiviruses infection
The human osteosarcoma cell line U-2 OS cells were purchased from Shanghai YiYan Biotech Inc. (Shanghai, China), and cultured in HyClone-Dulbecco’s modified eagle medium (DMEM) which contained 10% fetal bovine serum at 37℃ in a humidified incubator with 5% CO2. To build TIMP2 and BAX stable over-expression cells, we designed and purchased Lenti-TIMP2-EGFP and Lenti-BAX-EGFP from GeneChem Company (Shanghai, China). Based on the standard procedure, the U-2 OS cells were seeded in six-well plates supplied with DMEM and 10% fetal bovine serum. The cells were infected with Lenti-TIMP2-EGFP or Lenti-BAX-EGFP and Lenti-vector control, and stable cell clones were selected in the presence of puromycin for 1 to 3 weeks. The gene and protein expressions of TIMP2 and BAX in the cell line were evaluated using real-time PCR, agarose gel electrophoresis and Western blot by GeneChem Company (Shanghai, China).
Cell counting and growth curve assay
To evaluate the colony formation, the TIMP2 and BAX infected U-2 OS cells were firstly cultured to the logarithmic growth stage. Then cells were planted into four plates at a density of 50/100/150/200cells per plate and culture for 2 weeks. Finally, the cells with crystal violet stain were counted using microscope and quantified using ImageJ (NIH, USA).
The cells growth curve was evaluated by using real-time cell analyzer (RTCA) xCELLigence system (ACEA Bioscience, San Diego, CA). According to the manufacturer’s instructions, the cells were transferred to RTCA E-plates with a density of 1 × 104 cells per well and incubated for 1 hour. Then the baseline was recorded, and the electrical impedance was measured for four days in each well. The cell index was calculated based on the electrical impedance which can reflect cell viability.