Identifying candidate genes that are critical for the survival of osteosarcoma cells
Project Achilles knocks out each gene individually by using a genome-scale CRISPR-Cas9 tool and identifying candidate genes that are critical for the survival of tumors. Taking advantage of RNA-Seq data and Project Achilles from osteosarcoma patients, we can pinpoint essential genes responsible for osteosarcoma malignancy. Candidate genes were defined as essential genes responsible for osteosarcoma malignancy with a CERES score of < − 1 across 75% of osteosarcoma cell lines (n = 8), including Saos-2, G-292, U-2 OS, SJSA-1, 143B, HuO9, OS252, and HS-Os-1. In total, 785 genes were finally found to be crucial for maintaining survival in the 8 osteosarcoma cell lines (Additional file 1: S1).
To identify which candidate genes among these 785 genes were dysregulated expression in osteosarcoma tissues, DEG (Differentially Expressed Genes) analyses were performed to compare osteosarcoma tissues (n = 44) with normal tissues (n = 5) in the GSE19276 database. Fifty-nine of 785 essential genes were significantly changed in osteosarcoma tissues (Fig. 1, Additional file 2: S2). Among these 59 genes, 14 genes were significantly down-regulated, and 45 genes were significantly up-regulated in osteosarcoma tissues.
To identify the biological roles of these essential genes in osteosarcoma, the enrichment with KEGG-pathway and GO-BP (biological process) was analyzed. From the results of KEGG-pathway enrichment, we found that these 59 genes are mainly enriched in the proteasome, RNA transport, DNA replication, Nucleotide excision repair, aminoacyl tRNA biosynthesis, ribosome biogenesis in eukaryotes, ubiquitin mediated proteolysis, mismatch repair, base excision repair and homologous recombination (Fig. 2A). In addition, in the enrichment results of GO-BP, these genes were mainly related to the cell cycle, ribonucleoprotein complex biogenesis, ribosome biogenesis, RNA localization, transcription coupled nucleotide excision repair, nucleotide excision repair, anaphase promoting complex depend catabolic process, regulation of hematopoietic progenitor cell differentiation, and positive regulation of RNA polymerase II transcripition preinitiation complex assembly (Fig. 2B). In additional, enrichment analysis was conducted for the 45 upregulated genes and the 14 downregulated genes, respectively. However, the number of down-regulated genes is too small to perform a reliable KEGG and GO analysis. The KEGG and GO analysis were performed on only upregulated genes. From the results of KEGG-pathway enrichment, we found that these 45 upregulated genes are also mainly enriched in the proteasome, RNA transport, etc., (Fig. 2C). In the enrichment results of GO-BP, these 45 upregulated genes were also mainly related to the cell cycle, ribonucleoprotein complex biogenesis, etc. (Fig. 2D). The results showed that these essential genes were mainly associated with the cell cycle-related signaling pathway in the osteosarcoma tissues.
Next, a protein–protein interaction (PPI) network of these 59 essential genes was generated using Metascape (Fig. 2E). After the analysis of molecular complex detection, three MCODE components were obtained from the PPI network (Fig. 2F). In addition, the enrichment analysis of every MCODE component is depicted in Fig. 2D, in which the top three enrichment results are displayed. Among the three MCODE components, the most significant was MCODE1, which included six genes: LARS (Leucyl-TRNA Synthetase 1), BYSL (Bystin Like), UBA1 (Ubiquitin Like Modifier Activating Enzyme 1), PSMC6 (Proteasome 26S Subunit, ATPase 6), AARS1 (Alanyl-TRNA Synthetase 1) and ANAPC10 (Anaphase Promoting Complex Subunit 10).
Construction of essential genes-based prognostic predictors
The prognosis of these 59 essential genes for osteosarcoma patients was further analyzed by using the univariate Cox regression analysis based on GSE21257 dataset. The results revealed that high expression of LARS was significantly negatively correlated with the prognosis of osteosarcoma patients, and low expression of DNAJC17 (DnaJ Heat Shock Protein Family (Hsp40) Member C17) was significantly negatively correlated with the prognosis of osteosarcoma patients, Next, we used the LASSO analysis of the tenfold cross-validation to analyze LARS and DNAJC17 genes (Fig. 3B), and constructed the risk score model that included LARS and DNAJC17 genes (Fig. 3C). The risk score was calculated as follows: Risk score = (LARS × 1.7203 − DNAJC17 × 1.7652). Each osteosarcoma patient was divided into low-risk and high-risk groups according to the risk score. The distribution of risk scores, survival status, and the mRNA expression of two genes in osteosarcoma patients were shown in Fig. 3D. Kaplan–Meier survival revealed that the low-risk group was correlated with a better prognosis for osteosarcoma patients (Fig. 3E). The result of the ROC analysis demonstrated that the risk score had a powerful ability to predict prognosis in osteosarcoma patients (1-year (AUC = 0.88), 3-year (AUC = 0.80), 5-year (AUC = 0.75); Fig. 3F).
To further verify the prognostic model, we downloaded 37 osteosarcoma samples with complete clinical information from GSE39055 database. Each patient was brought into the previous prognostic model to calculate the risk score. Patients were divided into high-risk and low-risk groups. Kaplan–Meier curve analysis showed that osteosarcoma patients with low-risk scores had a better overall survival than those in the high-risk-score group (Fig. 4A). The result of the ROC analysis demonstrated that the risk score had a powerful ability to predict prognosis in osteosarcoma patients (1-year (AUC = 0.80), 3-year (AUC = 0.74), 5-year (AUC = 0.70); Fig. 4B).
Next, a nomogram model-independent predictors of OS, including risk score, gender, age, and tumor grade, was developed. The nomogram model provides a visual statistical predictive tool for the survival of osteosarcoma patients, allowing the risk-scoring model to be applied to clinical practice (Fig. 5A). The C‐index for the nomogram is 0.75, indicating the nomograph we built is reliable and accurate (Fig. 5B).
Analysis of drug sensitivity associated with LARS and DNAJC17
The GSCA database was exploited to analyze the relationship between the drug sensitivity and the expression of LARS and DNAJC17 based on the data from the Cancer Drug Sensitivity Genomics Database (GDSC) (Fig. 6, Additional file 1: S3). The result revealed that high LARS expression was associated with higher drug sensitivity to BRD-K30748066, Tozasertib, BI-2536, GSK461364, BRD-K70511574, Triazolothiadiazine, Rigosertib, FQI-1, KX2-391 and SB-225002 (R < − 0.3, p < 0.05). There was no significant correlation between the expression of DNAJC17 and the drug sensitivity.
Knockdown of LARS inhibits osteosarcoma cell proliferation
LARS and DNAJC17 are highly expressed in osteosarcoma. High expression of LARS is an unfavorable prognostic factor in patients with osteosarcoma. However, high expression of DNAJC17 is a favorable prognostic factor in patients with osteosarcoma. Consequently, we were interested in the influence of LARS on osteosarcoma. As noted earlier in this article, a CERES score < − 1 was defined as an essential gene for tumor cell survival. The CERES score of LARS in eight osteosarcoma cell lines was shown in Fig. 7A. The result indicated that the CERES score of LARS was all less than -1 in the above cells. Next, GSEA was performed to identify the differentially activated signaling pathways in the high LARS expression. The GSEA plot showed that high LARS expression positively correlated with the cell cycle (Fig. 7B). These results suggested that LARS may affect the proliferation of osteosarcoma cells. Then, LARS-targeting shRNAs (LARS-shRNA1 and LARS-shRNA2) were transfected into Saos-2 and U-2 OS cells, and LARS was successfully knocked down (Fig. 7C–F). CCK-8 analysis showed that LARS knockdown strongly inhibited the proliferation of Saos-2 and U-2 OS cell lines (Fig. 7G and H).
Evaluation of associations among LARS and immune cells
Our GSEA results showed that high expression of LARS was associated with down-regulation of the B cell receptor signaling pathway and T cell receptor signaling pathway (Fig. 8A). To further validate LARS as a potential immune influencer, we calculated the degree of immune infiltration through MCP-counter to estimate the tendency of residential immune cells to vary with expression changes in LARS. The results indicated that LARS expression is significantly negatively correlated with the infiltration levels of Monocytic lineage (R = − 0.42, p < 0.01), T cells (R = − 0.35, p < 0.05), and Fibroblasts (R = − 0.31, p < 0.05) in osteosarcoma (Fig. 8B). In addition, osteosarcoma patients were classified into low-group and high-group LARS subtypes based on the median LARS expression level. ESTIMATE results indicated that patients with lower LARS expression had a significantly higher stromal score and higher immune score, but with lower tumor purity relative to patients in high-group LARS (Fig. 8C and D).
The level of LARS protein was highly expressed in osteosarcoma
Finally, the protein expression of LARS in osteosarcoma was detected via IHC staining analysis. LARS protein expression was demonstrated significantly increased in 18 osteosarcoma tissues compared to the paired adjacent normal tissues (Fig. 9A–C).