Identification of MAP3K15 as a potential prognostic biomarker and correlation with immune infiltrates in osteosarcoma

Introduction

Osteosarcoma (OS) is a type of primary malignant tumor that commonly occurs in children and adolescents (1). Although its incidence is not high, OS easily relapses and metastasizes, and has a 5-year survival rate of approximately 65% (1,2). Currently, the treatment for OS mainly consists of 3 parts: preoperative chemotherapy, surgical resection, and postoperative radiotherapy (3,4). Surgical resection is still the primary treatment for OS (3). However, traditional medicines are unsatisfactory in reducing the recurrence and metastasis of OS (5,6).

Recently, immunotherapy has made breakthroughs in the field of cancer therapy. Immune checkpoint inhibitors, such as PD-1, PD-L1, and CTAL4, have shown significant antitumor effects on various cancers (7-9). Additionally, these immune checkpoint inhibitors have been shown to play substantial roles in the treatment of OS. Despite this, immunotherapy is still not beneficial to many patients because of target mismatch and drug resistance (10,11). Thus, there is an urgent need to identify novel immune-related biomarkers to better guide clinical immunotherapy for OS.

The tumor microenvironment (TME) is the cellular environment that includes endothelial cells, mesenchymal cells, immune cells, inflammatory mediators, and extracellular matrix molecules (12-14). The tumor can influence the microenvironment by releasing extracellular signals, promoting tumor angiogenesis and inducing peripheral immune tolerance, while the immune cells in the microenvironment can affect the growth and evolution of cancerous cells. the study of predictive biomarkers, which may fully address the complexities of the biology, will promote the development of therapies tailored to individual patients. Thus, it is important to understand the role of different components of the TME in the treatment and prevention of osteosarcoma. Further studies on the complexity and systemic characteristics of the TME are greatly needed to identify useful biomarkers. Based on the ESTIMATE and CIBERSORT algorithms, gene expression matrix data can be converted into an immune score and proportions of infiltrating immune cells, which can be used to study the immune activity of the TME (15-17). With the help of these 2 algorithms, many studies have explored immunity in the TME of lung cancer (18), breast cancer (19), and liver cancer (20). These studies provide direction for tumor immunotherapy by exploring the immune-related genes and immune cells that play crucial roles in the TME. Investigating the TME using these algorithms may provide new insights into immunotherapy for OS.

In this study, we collected the expression profiles and clinical information of OS patients using public databases, namely the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database and the GSE39058 dataset from the Gene Expression Omnibus (GEO) database. It was found that the immune score of the TME is statistically significant for the prognosis of patients with OS. Subsequently, we developed a prognostic signature based on the multivariate Cox regression and identified mitogen-activated protein kinase kinase kinase 15 (MAP3K15) as an important prognostic gene. We estimated infiltrating immune cells by the CIBERSORT algorithm and investigated the relationship between the prognostic signature and tumor immune infiltrating cells. Ultimately, a positive correlation was found between MAP3K15 expression and several infiltrating immune cells, which is expected to provide novel ideas for OS immunotherapy. In addition, OS patients with higher MAP3K15 expression had a poorer prognosis. We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/atm-21-3181).

Methods

Samples

A total of 48 OS tissues were obtained from patients at Sun Yat-sen Memorial Hospital of Sun Yat-sen University and The Eighth Affiliated Hospital of Sun Yat-sen University between December 2010 and March 2018. None of the patients underwent chemotherapy or radiotherapy before surgery. The tissues were immediately snap-frozen and stored at −80 °C in an ultralow temperature freezer until further analysis. The study was approved by the Ethics Committee of Sun Yat-sen Memorial Hospital of Sun Yat-sen University and The Eighth Affiliated Hospital of Sun Yat-sen University (Approval no.: 2019r013; Approval date: 1/3/2019). Informed consent forms were signed by all individuals prior to acquiring and studying their tissues. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Obtaining expression data from publicly available data

OS primary data were downloaded from the TARGET database and GEO database. Accession number GSE39058 was chosen from GEO for data analysis, and it contained a total of 164 samples, of which 46 included mRNA profiling data. Additionally, clinical outcomes were included in the dataset. Then, the microarray data of mRNA expression were normalized via a log2 scale transformation, and the average expression was used when a gene was detected with more than 1 probe. In addition, the RNA-seq data and relevant clinical information were downloaded from the TARGET database (https://ocg.cancer.gov/programmes/target). Both fragments per kilobase of transcript per million mapped reads (FPKM) and counts were used for further data processing. This study followed the publication guidelines of the GEO and TARGET databases.

Identification of the TME and infiltrating immune cells

First, the ESTIMATE algorithm (15) was performed to calculate the immune score and stromal score. According to the median values of the stromal/immune scores, we divided OS patients into a high-risk group and a low-risk group. Then, Kaplan-Meier survival curves were used to analyze the correlation between the clinical outcomes and the stromal/immune scores.

Additionally, the characterization of the cellular composition of OS was analyzed based on the CIBERSORT algorithm (17). A heatmap was used to visualize the proportions of immune cells in the low- and high-risk immune score groups. In addition, we further applied Spearman’s correlation coefficient to verify the relationship between immune/stromal scores and immune cells. Both the ESTIMATE and CIBERSORT algorithms were performed based on R 4.0.2 and RStudio (version 1.2.5001).

Differential expression analysis and functional analysis

The Limma R package (21) was applied to identify the differentially expressed immune genes (DEIGs) between the low and high immune score groups. Then, the common DEIGs were obtained from GSE39058 and TARGET-OS. Here, |log₂FC| ≥1 and P<0.05 were set as the cutoff values.

Furthermore, to explore the potential functions of these DEIGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed via the Metascape database (http://metascape.org). Additionally, overexpressed and underexpressed immune gene groups were input, and the protein-protein interaction (PPI) network and molecular complex detection (MCODE) were used to further identify the densely connected regions. P<0.05 was used as a cutoff value.

Construction of the prognostic immune gene signatures

TARGET-OS and GSE39058 were used to establish a prognostic immune gene signature, with TARGET-OS being the training set and GSE39058 being the validation set. Univariate Cox regression analysis was first performed to identify survival-related immune genes via the SangerBox database (http://sangerbox.com/). Additionally, a robust likelihood-based survival analysis was used to increase the reliability and accuracy, allowing for iterative forward selection of multiple sets of genes based on the partial likelihood of the Cox model.

Then, multivariate Cox regression was implemented to explore the independent prognostic predictor genes. Additionally, the survival receiver operating characteristic (ROC) curve was obtained from SangerBox, then the area under the ROC curve (AUC) was used to assess the efficiency of the predictive signature. P<0.05 was considered statistically significant.

Correlation analysis between immune cells and the prognostic signature

The proportions of 22 infiltrating immune cells were obtained according to the CIBERSORT algorithm. The Wilcoxon signed-rank test (2 groups) or Kruskal-Wallis rank sum test (more than 2 groups) were conducted to calculate the correlation between immune cells and clinical characteristics. Additionally, a log-rank test was used to evaluate the association between immune cells and the prognostic signature. P<0.05 was regarded as statistically significant.

Quantitative polymerase chain reaction (qPCR) assay

Extraction of total RNA was performed with TRIzol reagent (Life Technologies) following the manufacturer’s instructions, and the extracted RNA was then quantified. Next, qPCR was performed using a BeyoFast™ SYBR Green One-Step qRT-PCR Kit (Beyotime, D7268S). The expression of MAP3K15 was normalized to that of GAPDH. Relative expression levels were calculated with the 2^−∆∆Ct method.

Immunohistochemistry validation

Immunohistochemistry (IHC) was performed according to the antibody supplier instructions. The slices of the clinical samples were incubated with a primary antibody against MAP3K15 at different dilution ratios overnight at 4 °C. Images were captured at an appropriate magnification under a microscope (Nikon Microsystems, Shanghai, China). The antibody used in the study was anti-MAP3K15 (SAB, #45534).

Statistical Analysis

All statistical analyses were performed using R-version 4.0.3. The differential analysis was performed using the “limma R” package. The development of prediction signature was performed by “survival R” package, “survminer R” package, “caret R” package, “glmnet R” package, “time ROC R” package. The assessment of TME and infiltrated immune were calculated by the “ESTIMATE R” and “CIBERSORT R” package. The correlation plot was drawn by the “ggpubr R” package.

Results

Immune and stromal scores were evaluated by the ESTIMATE algorithm

First, 86 OS patients with gene expression profiles and clinical characteristics were identified in the TARGET-OS dataset. The ESTIMATE algorithm was then performed to calculate the immune and stromal scores of tumor tissues. Depending on the calculations, immune scores were distributed between −1,842.4 and 2,013.6, while stromal scores ranged from −2,124.1 to 3,901.0. To explore the prognostic significance of immune and stromal scores, we divided the OS patients into high and low immune score groups according to their scores. As shown in the Kaplan-Meier curves, although the survival differences in estimated group was not statistically significant (P=0.05), the stromal group (P=0.014) and immune group (P=0.0019) were positively correlated with prognosis, indicating that both were positive factors in the prognosis of OS patients (Figure 1A,B,C).

Figure 1 The prognostic significance of the estimated immune and stromal scores in the TARGET-OS and GSE39058 datasets. (A) Overall survival according to the estimated score in the TARGET-OS dataset. (B) Overall survival according to the immune score in the TARGET-OS dataset. (C) Overall survival according to the stromal score in the TARGET-OS dataset. (D) Overall survival according to the estimated score in the GSE39058 dataset. (E) Overall survival according to the immune score in the GSE39058 dataset. (F) Overall survival according to the stromal score in the GSE39058 dataset. TARGET, Therapeutically Applicable Research To Generate Effective Treatments; OS, osteosarcoma.

Additionally, the same processes were performed on the GSE39058 dataset. According to the ESTIMATE calculation, the immune scores ranged from −1,116.8 to 1,676.5, while stromal scores were between −2,677.4 and 2,916.059471. However, when the groups were divided into high and low immune score groups according to their scores, the results showed that there was no statistically significant difference in the prognosis between the immune group (P=0.23) and the stromal group (P=0.84) (Figure 1D,E,F). Although the survival analysis of immune score in GSE39058 was not significant, the curves of high and low immune score groups were divergent, which might be due to insufficient samples.

DEIGs and functional analysis

To explore the correlation between immune genes and immune scores, we compared 46 OS patients in the GSE39058 dataset and 86 patients in the TARGET-OS dataset. A total of 1,446 DEIGs were discovered in GSE39058, including 729 genes that were upregulated and 735 that were downregulated. In addition, a total of 565 DEIGs were identified in the TARGET-OS dataset, with 502 genes that were upregulated and 63 genes that were downregulated. Then, 100 commonly upregulated genes and 3 commonly downregulated genes were extracted via a Venn diagram (Figure 2A,B,C).

Figure 2 Overview of the DEIGs in OS. (A) A heatmap of the expression profiles in the TARGET-OS dataset. (B) A heatmap of the expression profiles in the GSE39058 dataset. (C) A Venn diagram showing the common DEIGs in both datasets. (D) Bar plots displaying GO and KEGG analysis results based on Metascape. (E) Network display of GO and KEGG analysis results based on Metascape. (F) MCODE results obtained from the PPI network. DEIG, differentially expressed immune gene; OS, osteosarcoma; TARGET, Therapeutically Applicable Research To Generate Effective Treatments; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction.

Subsequently, GO enrichment analysis indicated that the common DEIGs were enriched in various immune processes, such as leukocyte proliferation, regulation of leukocyte differentiation, negative regulation of leukocyte activation, regulation of innate immune response, and leukocyte-mediated cytotoxicity, while KEGG pathway analysis suggested that osteoclast differentiation, natural killer (NK) cell-mediated cytotoxicity, and the adaptive immune system were the most pivotal pathways affected by these genes (Figure 2D,E). Thus, these dysregulated immune genes might be the basis of the changes in the TME and might enhance the mechanism of immune escape.

Furthermore, the PPI network and MCODE plugin based on the Metascape database identified the significant modules in these DEIGs. Module 1 included 33 edges and 9 nodes, and involved MAP3K15, CDC247, and HLA-DRB1. Module 2 included 6 edges and 4 nodes, and involved CCL23, CCL5, GPR55, and LNP3 (Figure 2F).

Identification of a prognosis-related gene

Since 103 common DEIGs were obtained from the previous step, univariate Cox regression analysis was first performed to filter the survival-related genes. According to the median expression levels of immune genes, OS patients were divided into high and low expression groups. As shown in Figure 3, 3 and 12 genes were identified to be correlated with prognosis in the GSE39058 and TARGET-OS datasets, respectively.

Figure 3 Survival curves of the DEIGs were created using Kaplan-Meier curves. (A) DEIGs with prognostic value in the TARGET-OS dataset. (B) DEIGs with prognostic value in the GSE39058 dataset. DEIG, differentially expressed immune gene; TARGET, Therapeutically Applicable Research To Generate Effective Treatments; OS, osteosarcoma.

Then, we noticed that only MAP3K15 was common in the results from the 2 datasets. Thus, multivariate Cox regression was implemented to explore the independent prognostic value of MAP3K15. As shown in Figure 4, MAP3K15 was verified to be statistically significant. Moreover, ROC curves were used to determine the efficiency of the Cox regression, and the 5-year AUC was 0.65. In addition, the result was verified in GSE39058, with an AUC of 0.63 (Figure S1). Therefore, the expression of MAP3K15 could serve as an independent prognostic signature to predict OS patients’ clinical outcomes, and may provide novel insights into the changes in the TME and the mechanism of immune escape.

Figure 4 Kaplan-Meier plot of MAP3K15 and ROC curve of the prognostic immune signature. (A) The predictive immune model constructed from MAP3K15. (B) ROC curves of the immune signature with the corresponding AUCs indicated on the side. ROC, receiver operating characteristic; AUC, area under the receiver operating characteristic curve.

Immune landscapes evaluated by CIBERSORT

There are a considerable number of tumor cells and immune cells in the TME. Understanding the immune landscape of the tumor may reveal the underlying mechanism of refractory tumors. Thus, the CIBERSORT algorithm was applied to explore the proportions of immune cells in the GSE39058 and TARGET-OS datasets based on the gene expression data. Here, analysis of the GSE39058 and TARGET-OS datasets was processed separately. Although the proportions of immune cells suggested different immune landscapes in the 2 datasets, CD4⁺ naive T cells, CD4⁺ memory resting T cells, M0 macrophages, M2 macrophages, and naive B cells accounted for a large proportion of infiltrating immune cells in OS (Figure 5A,B). Moreover, violin plots were used to visualize the immune cell subsets between the high and low immune score groups (Figure 5C,D), and the results showed that the proportions of memory B cells, naive T cells, M2 macrophages, and activated NK cells were significantly increased.

Figure 5 Immune cell landscapes in different immune score groups. (A) The proportion of immune cells in the TARGET-OS dataset. (B) The proportion of immune cells in the GSE39058 dataset. (C) High immune score and low immune score groups in the TARGET-OS dataset. (D) High immune score and low immune score groups in the GSE39058 dataset. (E) Proportions of memory B cells, CD4 T cells, activated NK cells, M0 macrophages, M1 macrophages, and M2 macrophages in different immune score groups (TARGET-OS). (F) Proportions of memory B cells, CD4 T cells, activated NK cells, M0 macrophages, M1 macrophages, and M2 macrophages in different immune score groups (GSE39058). *P<0.05; **P<0.01. TARGET, Therapeutically Applicable Research To Generate Effective Treatments; OS, osteosarcoma.

Additionally, to investigate the rationality of our predictive signature, we further calculated the proportions of immune cells in the high- and low-risk immune score groups (Figure 5E,F). The results demonstrated that the proportions of CD4⁺ T cells, M1 macrophages, and M2 macrophages were significantly increased. Thus, the differences between the high- and low-risk score groups may provide insights into the mechanism of immune microenvironment regulation.

Correlations between the prognostic signature and infiltrating immune cells

Since the results above revealed that MAP3K15 might play important roles in the TME, we questioned the association between MAP3K15 and the relevant immune cells. According to the calculation of immune cells by the CIBERSORT algorithm, a correlation between MAP3K15 and infiltrating immune cells was found. As shown in Figure 6A,B,C, the expression of MAP3K15 and CD4⁺ T cells, M0 macrophages, and M1 macrophages were positively correlated, which indicated that MAP3K15 in involved in immune cytotoxicity.

Figure 6 Correlations between the prognostic signature and infiltrating immune cells. (A,B,C) The expression of MAP3K15 and CD4 T cells, M0 macrophages, and M1 macrophages were positively correlated. (D) MAP3K15 was positively correlated with the immune score.

Furthermore, according to the ESTIMATE algorithm, we found that MAP3K15 was positively correlated with the immune score, which further suggested the significance of MAP3K15 in the TME (Figure 6D).

Validation

To verify the results obtained by bioinformatics analysis, we performed a qPCR assay of MAP3K15 expression in OS tissues from 48 patients. Primer sequences are shown in Table 1. The patients were divided into 2 groups, the high expression and low expression groups, according to MAP3K15 expression. Results showed that patients with high MAP3K15 expression had a significantly poorer prognosis (Figure 7A). Moreover, the IHC results further confirmed the expression trend of MAP3K15 protein in these samples (Figure 7B).

Table 1 Primer sequence references for MAP3K15 and GAPDH
Full table

Figure 7 MAP3K15 expression validation and prognosis analysis. (A) The OS rates of patients with high and low expression of MAP3K15 were evaluated by the log-rank method. (B) Representative histopathological images of MAP3K15 expression in OS. Brown staining of cell indicates positive staining and the lack of brown staining indicates negative staining. Scale bar =100 µm for ×200, scale bar =50 µm for ×400. OS, osteosarcoma.

Discussion

Recently, tumor immunotherapy has attracted an increasing amount of attention, and it is considered one of the most likely treatment methods to improve the survival rate of cancer patients. Immunotherapy has been successful in a variety of malignant tumors, including lung cancer, and several immunotherapy drugs have been approved by the Food and Drug Administration (FDA) (22,23). However, the application of immunotherapy in OS has been slow due to the complexity of the immune system and the subtle differences in the tumor-specific microenvironment, as well as the fact that the tumor itself may resist immunotherapy in various ways (24). It is necessary to determine predictive biomarkers with higher sensitivity and specificity to stratify patients and to maximize the benefits for patients. In the TME, the interactions among tumor cells, immune cells, and stromal cells promote tumor proliferation and metastasis (25,26). Therefore, a clear understanding of the complex TME and the complex immune characteristics of tumors and immune cells may be helpful for guiding clinical immunotherapy.

The ESTIMATE and CIBERSORT algorithms are currently widely used to study the infiltration of immune cells and the proportions of cell subtypes in the TME. These 2 algorithms can transform the information of the gene expression matrix into an immune score and immune cell infiltration ratio. For instance, Lai et al. (27) studied the TME in renal cell carcinoma and identified effective biomarkers. Moreover, other tumor-related studies, such as lung cancer (28,29) and liver cancer (30,31), were also reported to use these algorithms. Interestingly, Hong et al. (32) studied the TME of OS using bioinformational methods, but no validation was performed in OS tissue. Therefore, further studies to understand the immune microenvironment of OS will be helpful to reveal the mechanism of tumor occurrence and discover effective biomarkers.

We searched the currently available public databases for OS, and the TARGET-OS and GSE39058 datasets were used in this study. Through the ESTIMATE algorithm, TARGET-OS showed that the immune score of the immune microenvironment of OS was closely related to prognosis. That is, the higher the immune score, the better the prognosis. Although no statistically significant differences were found in GSE39058, it can be concluded from the survival curve that there was a significant difference between the high and low immune score groups. Therefore, we speculated that this might be caused by the insufficient number of patients.

Then, we obtained the DEIGs by analyzing the difference between the high and low immune score groups and determined the intersection of DEIGs in the 2 datasets. The GO results revealed that these DEIGs were involved in various key immune responses, such as leukocyte proliferation, regulation of leukocyte differentiation, negative regulation of leukocyte activation, regulation of innate immune response, and leukocyte-mediated cytotoxicity, while KEGG pathway analysis suggested that osteoclast differentiation, NK cell-mediated cytotoxicity, and the adaptive immune system were the most pivotal pathways affected by these genes. These findings provide a direction for us to further study the mechanism of the immune response. Using univariate and multivariate regression analyses, we found that only MAP3K15 was statistically significant in both databases. MAP3K15, also known as ASK3, is a member of the mitogen-activated protein kinase (MAPK) family. To our knowledge, this is the first report of MAP3K15 in tumor. However, some studies have found that the MAPK family members function in a protein kinase signal transduction cascade, where an activated MAPK kinase kinase (MAP3K) phosphorylates and activates a specific MAPK kinase (MAP2K), which then activates a specific MAPK. This MAP3K protein plays an essential role in apoptotic cell death triggered by cellular stresses (33). This result suggests that MAP3K15 may play a significant role in tumor initiation and progression.

Subsequently, the CIBERSORT algorithm further identified the proportions of infiltrating immune cells in the TME of OS. We found that memory B cells, naive T cells, M2 macrophages, and NK cell-activated cells were the main infiltrating cells. Then, we discussed the correlations between these immune cells and MAP3K15, as well as the correlations between these immune cells and the high- and low-risk groups. Finally, we found that MAP3K15 can predict the prognosis of patients with OS, and MAP3K15 is closely related to CD4⁺ T cells, M0 macrophages, and M1 macrophages. These results may provide new ideas for immunotherapy in the clinic.

We noticed that Hong et al. (32) also studied the microenvironment of OS. Although we both studied the TME of OS, we believe that the 2 studies are different. In this research, we found a common DEIG (MAP3K15) by intersecting the analysis results, while in the study of Hong et al., they determined that SIGLEC7 and SP140 were DEIGs. This is reasonable as different research strategies were utilized. Moreover, we performed Cox regression analysis and analyzed the correlation between MAP3K15 and immune cells. Specifically, no validation assays were performed in Hong et al.’s study, and they did not further study the relationship between genes and immune cells.

Immune cell infiltration and the TME are hot research topics at present, and the development of bioinformatics methods has further improved our understanding of tumor molecular biology. However, there are still many limitations in this study. Although we validated the expression of MAP3K15 using qPCR and IHC and found that patients with high MAP3K15 expression had a significantly poorer prognosis, the correlation between infiltrating immune cells and MAP3K15 expression was not tested because of restrictions associated with our experimental conditions. Since we found a potential association between MAP3K15 expression and immunity, further mechanisms need to be explored.

Taken together, our study found that MAP3K15, whose expression level is closely related to immune activity in tumors, is a critical immune-related biomarker, and our findings may provide a basis for OS immunotherapy.

Acknowledgments

Funding: This study was financially supported by the National Natural Science Foundation of China (81871750, 81971518), the Key Realm R&D Program of Guangdong Province (2019B020236001), and the Shenzhen Key Medical Discipline Construction Fund (ZDSYS20190902092851024).

Footnote

Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://dx.doi.org/10.21037/atm-21-3181

Data Sharing Statement: Available at https://dx.doi.org/10.21037/atm-21-3181

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/atm-21-3181). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was approved by the Ethics Committee of Sun Yat-sen Memorial Hospital of Sun Yat-sen University and The Eighth Affiliated Hospital of Sun Yat-sen University (Approval No.: 2019r013; Approval date: 1/3/2019). Informed consent forms were signed by all individuals prior to acquiring and studying their tissues. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin 2018;68:7-30. [Crossref] [PubMed]
2. Pingping B, Yuhong Z, Weiqi L, et al. Incidence and mortality of sarcomas in Shanghai, China, during 2002-2014. Front Oncol 2019;9:662. [Crossref] [PubMed]
3. Ritter J, Bielack SS. Osteosarcoma. Ann Oncol 2010;21:vii320-5. [Crossref] [PubMed]
4. Yang S, Ye Z, Wang Z, et al. High mobility group box 2 modulates the progression of osteosarcoma and is related with poor prognosis. Ann Transl Med 2020;8:1082. [Crossref] [PubMed]
5. Bishop MW, Janeway KA, Gorlick R. Future directions in the treatment of osteosarcoma. Curr Opin Pediatr 2016;28:26-33. [Crossref] [PubMed]
6. Dass CR, Choong PF. Gene therapy for osteosarcoma: steps towards clinical studies. J Pharm Pharmacol 2008;60:405-13. [Crossref] [PubMed]
7. Yoshida K, Okamoto M, Aoki K, et al. A review of T-cell related therapy for osteosarcoma. Int J Mol Sci 2020;21:4877. [Crossref] [PubMed]
8. Yoshida K, Okamoto M, Sasaki J, et al. Anti-PD-1 antibody decreases tumour-infiltrating regulatory T cells. BMC Cancer 2020;20:25. [Crossref] [PubMed]
9. Thanindratarn P, Dean DC, Nelson SD, et al. Advances in immune checkpoint inhibitors for bone sarcoma therapy. J Bone Oncol 2019;15:100221 [Crossref] [PubMed]
10. Le Cesne A, Marec-Berard P, Blay JY, et al. Programmed cell death 1 (PD-1) targeting in patients with advanced osteosarcomas: results from the PEMBROSARC study. Eur J Cancer 2019;119:151-7. [Crossref] [PubMed]
11. Zheng B, Ren T, Huang Y, et al. PD-1 axis expression in musculoskeletal tumors and antitumor effect of nivolumab in osteosarcoma model of humanized mouse. J Hematol Oncol 2018;11:16. [Crossref] [PubMed]
12. Wu T, Dai Y. Tumor microenvironment and therapeutic response. Cancer Lett 2017;387:61-8. [Crossref] [PubMed]
13. Pan H, Yu T, Sun L, et al. LncRNA FENDRR-mediated tumor suppression and tumor-immune microenvironment changes in non-small cell lung cancer. Transl Cancer Res 2020;9:3946-59. [Crossref]
14. Cersosimo F, Lonardi S, Bernardini G, et al. Tumor-associated macrophages in osteosarcoma: from mechanisms to therapy. Int J Mol Sci 2020;21:5207. [Crossref] [PubMed]
15. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
16. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
17. Gentles AJ, Newman AM, Liu CL, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med 2015;21:938-45. [Crossref] [PubMed]
18. Mao S, Li Y, Lu Z, et al. Systematic profiling of immune signatures identifies prognostic predictors in lung adenocarcinoma. Cell Oncol (Dordr) 2020;43:681-94. [Crossref] [PubMed]
19. Deng L, Lu D, Bai Y, et al. Immune profiles of tumor microenvironment and clinical prognosis among women with triple-negative breast cancer. Cancer Epidemiol Biomarkers Prev 2019;28:1977-85. [Crossref] [PubMed]
20. Liu F, Qin L, Liao Z, et al. Microenvironment characterization and multi-omics signatures related to prognosis and immunotherapy response of hepatocellular carcinoma. Exp Hematol Oncol 2020;9:10. [Crossref] [PubMed]
21. Law CW, Alhamdoosh M, Su S, et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Res 2016;5:ISCB Comm J-1408.
22. Subbiah V, Solit DB, Chan TA, et al. The FDA approval of pembrolizumab for adult and pediatric patients with tumor mutational burden (TMB) ≥10: a decision centered on empowering patients and their physicians. Ann Oncol 2020;31:1115-8. [Crossref] [PubMed]
23. Wu Y, Mou H. Complete remission in a patient with advanced renal pelvis carcinoma with lung metastasis treated with durvalumab immunotherapy: a case report. Ann Palliat Med 2020;9:2408-13. [Crossref] [PubMed]
24. Pratt HG, Justin EM, Lindsey BA. Applying osteosarcoma immunology to understand disease progression and assess immunotherapeutic response. Adv Exp Med Biol 2020;1258:91-109. [Crossref] [PubMed]
25. Cao R, Yuan L, Ma B, et al. Tumour microenvironment (TME) characterization identified prognosis and immunotherapy response in muscle-invasive bladder cancer (MIBC). Cancer Immunol Immunother 2021;70:1-18. [Crossref] [PubMed]
26. Giraldo NA, Sanchez-Salas R, Peske JD, et al. The clinical role of the TME in solid cancer. Br J Cancer 2019;120:45-53. [Crossref] [PubMed]
27. Lai Y, Tang F, Huang Y, et al. The tumour microenvironment and metabolism in renal cell carcinoma targeted or immune therapy. J Cell Physiol 2021;236:1616-27. [Crossref] [PubMed]
28. Zhang L, Chen J, Cheng T, et al. Identification of the key genes and characterizations of Tumor Immune Microenvironment in Lung Adenocarcinoma (LUAD) and Lung Squamous Cell Carcinoma (LUSC). J Cancer 2020;11:4965-79. [Crossref] [PubMed]
29. Liu Z, Wan Y, Qiu Y, et al. Development and validation of a novel immune-related prognostic model in lung squamous cell carcinoma. Int J Med Sci 2020;17:1393-405. [Crossref] [PubMed]
30. Liu X, Niu X, Qiu Z. A five-gene signature based on stromal/immune scores in the tumor microenvironment and its clinical implications for liver cancer. DNA Cell Biol 2020;39:1621-38. [Crossref] [PubMed]
31. Zhang FP, Huang YP, Luo WX, et al. Construction of a risk score prognosis model based on hepatocellular carcinoma microenvironment. World J Gastroenterol 2020;26:134-53. [Crossref] [PubMed]
32. Hong W, Yuan H, Gu Y, et al. Immune-related prognosis biomarkers associated with osteosarcoma microenvironment. Cancer Cell Int 2020;20:83. [Crossref] [PubMed]
33. Naguro I, Umeda T, Kobayashi Y, et al. ASK3 responds to osmotic stress and regulates blood pressure by suppressing WNK1-SPAK/OSR1 signaling in the kidney. Nat Commun 2012;3:1285. [Crossref] [PubMed]

(English Language Editor: C. Betlzar)