All the data in this study come from the target database (https://ocg.cancer.gov/programs/target), and the GEO database (https://www.ncbi.nlm.nih.gov/geo/), the target database contains the mRNA expression profile and clinical data of 98 osteosarcoma patients were obtained from the GEO database GSE16088 and GSE16091 data sets. The mRNA expression profile data and clinical data of 6 normal controls and 48 osteosarcoma patients were obtained. Normalize and debatch the obtained GEO data, standardize the samples, and view the clustering of the sample groups through the PCA graph and the UMAP graph.

The DESeq2 package of R software was used for expression profile difference analysis as well as gene volcano and heat map plotting. Corrected P values and |log2FC| were calculated, with P<0.05 and |log2FC|≥1.0 for differentially expressed genes. The Gene Ontology (GO) analysis was visualized using the ClusterProfiler package in R software, and a false discovery rate <0.01 was statistically significant.

The GSE16091 dataset included expression profiles, and survival events of 34 osteosarcoma tissues, which were suitable for the construction of weighted gene co-expression networks. GSE160306 and GSE160308 gene expression data matrices were constructed using the R package “ WGCNA”, and the top 25% of genes with the largest differences in the samples were selected as the input dataset for the subsequent WGCNA. A sample hierarchical clustering approach was used to detect and remove abnormal samples before selecting the appropriate soft threshold function. The neighbor-joining matrix and topological overlap matrix were constructed to calculate the corresponding dissimilarities and complete the identification of gene trees and modules. The minimum module size was 20. then the highly similar modules were merged by clustering and fusing module feature genes. The degree of difference is less than 0.25. Average linkage hierarchical clustering was established, similar genes were divided into a module, and osteosarcoma grading was used as phenotype data to screen the optimal module.

The impact of prognostic survival of osteosarcoma was assessed by overall survival (OS). All survival analyses were performed using COX multiple survival regression model analysis as the primary method and the Kaplan-Meier method as an adjunct. The endpoint of overall survival was defined as the time from the time of randomization to death from any cause, and all methods and R packages were performed using R software version v4.0.3. p < 0.05 was statistically significant.

The collection between gene expression and OS was assessed using LASSO COX analysis. A prognostic risk prediction model for osteosarcoma was built. Patients with osteosarcoma were divided into high and low-risk groups based on the median risk score. KM curves were plotted to compare OS between the high-risk and low-risk parts, and ROC survival analysis was performed using the R package SURVIVAL, and the “rmda” package was used for decision curve analysis.

The derived human osteoblast cell line hfob 1.19 and osteosarcoma cell line saos-2 were cultured in a 1:1 mixture of Ham’s F12 medium and Dulbecco’s modified Eagle’s medium containing 2.5 mM L -gu Aminoamide (ATCC, LGS standard, PL), supplemented with 100 U/mL penicillin, 100 U/mL streptomycin (Sigma-Aldrich, Warsaw, Poland), 0.3 mg/mL G418 (Sigma-Aldrich, Warsaw, Poland) and 10% fetal bovine serum (v/v, FBS, Gibco, Thermo Fisher Scientific, Warsaw, Poland). Cells were grown at 34°C, 5% CO₂. The pc DNA3.1 transfection plasmid was purchased from Genechem Co., Ltd. (Shanghai, China). Osaos-2 cells were seeded in 6-well plates (2×10⁵/well), and after incubation for 24 hours, they were transfected with LiPofectamine 2000 (corresponding plasmid 100 nM), and the operation was performed according to the instructions. Subsequent experiments were performed 48 hours after transfection.

Total RNA was extracted from each group of cells using TRIzol reagent (Invitrogen). cDNA was synthesized from RNA using a Primescript RT kit (Takara, Japan). The cDNA was amplified and quantified using SYBR Green mix (Finnish Finnzymes) in an Applied Biosystems 7500 instrument. The primer sequences used are shown in Table 1 . -2^ΔΔCt method could determine the gene expression.

Cells were collected, lysed in lysate (phosphatase inhibitor, protease inhibitor and PMSF), and protein quantification was performed by BCA method (Thermo Fisher Scientific, China). 10-20 μg protein was loaded on 8%-12% 30% acrylamide-Bis gel, then transferred through 0.22 μm pore PVDF membrane (Merck Millipore, USA), and blocked with 5% nonfat milk powder for 1 h., respectively, add primary antibodies (Rabbit polyclonal to IFIT1, 1:1000, ab236256; Rabbit monoclonal [EPR10009] to PARM-1, 1:1000, ab168369, abcam) and incubate overnight at 4°C. The next day, secondary antibodies (Goat Anti-Rabbit actin, 1:2000, ab8226, abcam) were added and incubated at room temperature for 1 h, and then luminescent solution (Thermo Fisher Scientific, China) was added for exposure and color development. Image J software was used for analysis, and the relative protein content was expressed as the gray value of the corresponding protein band/the gray value of the actin protein band, which was repeated three times.

Digest the cells and adjust to 1×106 cells/mL in culture medium without serum. Matrigel (BD Biosciences, San Jose, CA, USA) was added to the upper chamber and incubated at 37°C for 2 h. Add 100 µL of cell suspension to the upper chamber and 600 µL of cell culture medium with FBS to the bottom chamber. Subsequently, the cells were incubated at 37°Cand 5% CO₂ in an ambient environment for 24 hours. Detach the upper chamber, remove the cells, and fix with 4% paraformaldehyde. Cells were washed with phosphate-buffered saline (PBS) and stained with crystal violet. Stained cells were observed with a microscope (Olympus, Tokyo, Japan).

SPSS 24.0 and GraphPad Prism 8.0.1 software can process the data. Normally distributed measures are expressed as (x ± s), and comparisons between two parts were made by independent samples t-test, and P was a two-sided test, and differences were statistically significant at P < 0.05.

We downloaded the dataset of GSE16088 from the GEO database, which included 14 cases of osteosarcoma tissues and 6 cases of normal tissues, first, we analyzed the datasets, which showed good normalization and differences between samples ( Figures 1A–C ). We then analyzed the differentially expressed genes in normal and osteosarcoma tissues, and there were 5005 differentially expressed genes (2719 significantly up-regulated genes, 2286 down-regulated genes), Figures 1D, E show the differential gene volcano and heat map.

To identify genes associated with osteosarcoma survival orientation, we constructed co-expression networks by WGCNA in GSE16091, and in this study, we chose β=6 as the threshold to construct scaling-free networks ( Figures 2A, B ), and a total of 25 gene expression modules were identified after using merged dynamic tree cuts ( Figure 2C ), by constructing a random gene network map ( Figure 2D ). By calculating the correlation between module feature genes and clinical features, the skyblue3 module was most strongly correlated with the survival of osteosarcoma patients ( Figure 2E ), and the neighbor-joining matrix was clustered with the clinical features of the samples in Figure 2F , and the scatter plot of skyblue3 module genes in Figure 2G . We then obtained a total of 9 genes (FI27 ISG15, IFI6, IFI44L, PARM1, MX1, IFIH1, IFIT1, RSAD2) by taking the intersection of skyblue3 module genes and DEGs Figure 2H .

To explore the prognostic impact of intersecting genes of modular genes and DEGs on osteosarcoma, we did prognostic models based on LASSO analysis for these nine genes, obtained the value of the independent variable lambda and the coefficient of the independent variable ( Figure 3A ), and plotted the partial likelihood deviation versus log(λ) ( Figure 3B ). Patients with osteosarcoma were subsequently divided into high and low-risk groups ( Figure 3C ) by the hazard score Riskscore=(-0.8635)* IFIT1+(-0.481)* PARM1, and their survival status was presented in Figure 3C . Figure 3D shows the difference in survival between the high and low-risk parts (p=4.3e-3). Prognostic survival was predicted for 1, 3, and 5 years, and this showed good sensitivity and specificity ( Figure 3E ).

In order to further investigate the role of IFIT1 and PARM1 in patients with osteosarcoma, we analyzed the expression and prognostic effects of IFIT1 and PARM1 in patients with osteosarcoma, using the risk scores of IFIT1 and PARM1. IFIT1 expression was significantly increased in osteosarcoma tissues compared to normal tissues ( Figure 4A ), and patients with low IFIT1 expression had a better prognosis ( Figure 4B , p=0.07). PARM1 was lowly expressed in osteosarcoma tissues ( Figure 4C ), and low expression of IFIT1 predicted a poor prognosis for patients ( Figure 4D ).

To further verify the heterogeneity of IFIT1 and PARM1 in normal and osteosarcoma tissues, we cultured the osteoblast cell line hfob 1.19 and the osteosarcoma cell line saos-2, and detected the expression of IFIT1 and PARM1 by RT-qPCR, and the results showed that compared with osteoblasts, osteosarcoma cells had high expression of IFIT1 ( Figure 5A ), PARM1 was lowly expressed ( Figure 5B ). Then we detected the protein levels of IFIT1 and PARM1 in the osteoblast cell line hfob 1.19 and the osteosarcoma cell line saos-2 by WB. The results showed that the expression of IFIT1 was high and the expression of PARM1 was low in the osteosarcoma cell line saos-2 ( Figure 5C ). The above data show that IFIT1 is highly expressed and PARM1 is lowly expressed in osteosarcoma tissues. Subsequently, we overexpressed PARM1 in saos-2 cells, the transfection efficiency is shown in Figure 5D , and on this basis, we detected the cell migration ability. The results showed that overexpressing PARM1 in osteosarcoma cells significantly reduced the migration ability of cells ( Figure 5E , p<0.01).

In order to further verify the correctness of the constructed risk model, we constructed a prognostic risk model based on the target using 9 genes, FI27 ISG15, IFI6, IFI44L, PARM1, MX1, IFIH1, IFIT1, RSAD2, and obtained the independent variable lambda values ​​and the coefficients of the independent variable ( Figure 6A ), and plotted the partial likelihood deviation versus log(λ) ( Figure 6B ). Subsequently, patients with osteosarcoma were divided into high-risk group and low-risk group by the risk score Riskscore=(-0.0921)*PARM1+(-0.1418)*IFIH1 ( Figure 6C ), and their survival status was shown in Figure 3C . Figure 6D shows the difference in survival between high and low risk groups (P=0.00725). Prognostic survival was predicted for 1, 3, and 5 years, and the prognostic model showed good sensitivity and specificity ( Figure 6E ). Subsequently, we analyzed the effect of high or low expression of PARM1 on the prognosis of patients with osteosarcoma in the target database. Figure 7A shows the distribution of patients with osteosarcoma. The KM curve shows that patients with high PARM1 expression have a better prognosis ( Figure 7B ), and the ROC curve shows that patients with osteosarcoma have a better prognosis. Good sensitivity ( Figure 7C ).

In the previous data, we confirmed that low expression of PARM1 predicted poor prognosis in osteosarcoma patients. To further elaborate on the role of PARM1 in cancer, we performed a pan-cancer expression and prognosis analysis of PARM1. As shown in Figure 8A , PARM1 was lowly expressed in BLCA, BRCA, CESC, COAD, GBM, HNSC, KICH, KIRC, KIRP, LUAD, LUSC, PRAD, PEAD, THCA, and UCEC, and highly expressed in CHOL and LIHC. In addition, the expression level of PARM1 was closely correlated with the prognosis of COAD, KIRC, LUAD, and THYM, where low expression of PARM1 in COAD, KIRC, and LUAD predicted poor patient prognosis ( Figure 8B ).