Integrated Analysis of Multi-Omics Data to Establish a Hypoxia-Related Prognostic Model in Osteosarcoma

Abstract

Background:

Osteosarcoma (OS) is the most common malignant bone tumor in clinical practice, and currently, the ability to predict prognosis in the diagnosis of OS is limited. There is an urgent need to find new diagnostic methods and treatment strategies for OS.

Material and methods:

We downloaded the multi-omics data for OS from the TARGET database. Prognosis-associated methylation sites were used to identify clustered subtypes of OS, and OS was classified into 3 subtypes (C1, C2, C3). Survival analysis showed significant differences between the C3 subtype and the other subtypes. Subsequently, differentially expressed genes (DEGs) across subtypes were screened and subjected to pathway enrichment analysis.

Results:

A total of 249 DEGs were screened from C3 subtype to other subtypes. Metabolic pathway enrichment analysis showed that DEGs were significantly enriched to the hypoxic pathway. Based on univariate and multivariate COX regression analysis, 12 genes from the hypoxia pathway were further screened and used to construct hypoxia-related prognostic model (HRPM). External validation of the HRPM was performed on the GSE21257 dataset. Finally, differences in survival and immune infiltration between high and low risk score groups were compared.

Conclusion:

In summary, we proposed a hypoxia-associated risk model based on a 12-gene expression signature, which is potentially valuable for prognostic diagnosis of patients with OS.

Keywords

Osteosarcoma biomarkers hypoxia prognostic model methylation

Background

Osteosarcoma is a highly aggressive malignancy that is often fatal if not detected and treated early.¹ The 5-year survival rate of patients with osteosarcoma without distant metastases has now been increased to 60-70% percent in clinical practice using a mix of chemotherapy and surgical treatment options.² However, for patients with metastatic osteosarcoma, the 5-year survival rate is less than 30%.³ The capacity to predict prognosis at the time of osteosarcoma diagnosis is currently limited. At the moment of diagnosis, the presence or absence of metastases is the most important predictor of outcome.⁴ However, markers of metastasis are absolutely meaningless in patients with just localized disease. As a result, in order to achieve further progress in treatment and prognosis, more generalized predictive methods are needed.

Hypoxia is a characteristic of the great majority of solid tumors, in which tumor cell proliferation results in hypoxia due to an aberrant lack of blood supply to the tumor microvasculature.⁵ The absence of oxidative phosphorylation in a chronic hypoxic environment causes a drop in ATP, forcing cells to employ anaerobic respiration as a method of survival, resulting in an accumulation of lactate in the cells and a decrease in pH.⁶ Although extreme hypoxia causes tumor necrosis, moderate hypoxia close to the tumor center promotes tumor angiogenesis, cancer cell survival, and stem cell characteristics, encouraging tumor development, metastasis, and medication resistance.⁷ DNA methylation plays a significant role in defining tissue-specific gene expression and chromosomal instability through epigenetic regulation of gene expression.⁸ Overall hypomethylation plays a key function in the tumor cell hypoxia response system to adapt to hypoxia.⁹ Aberrant methylation has been connected to carcinogenesis, and in many tumor forms, variations in methylation patterns between malignancies and benign tissues have been observed.¹⁰ Several biomarkers related with DNA methylation have been revealed for early cancer diagnosis or prognostic prediction due to the highly conserved feature of DNA methylation.^11-14

In this study, we identified osteosarcoma methylation subtypes using prognostic methylation loci and further screened for differentially expressed genes (DEGs) between subtypes. Metabolic pathway enrichment analysis of the DEGs indicated that the hypoxic metabolic pathway was significantly enriched. Accordingly, we constructed a hypoxia-related prognostic model hypoxia-related prognostic model (HRPM) based on hypoxia pathway genes.

Material and Methods

Data acquisition

We downloaded gene expression, methylation and copy number variation (CNV) data for OS from the TARGET database (https://ocg.cancer.gov/programs/target). Clinical information of OS patients was also obtained from the TARGET database, and samples lacking survival information were excluded. Next, we matched only the sample labels shared between gene expression data, methylation data, CNV data and clinical data, and finally retained 83 OS samples. The GSE21257 dataset contains 53 samples of pre-chemotherapy osteosarcoma patients, we obtained their gene expression matrix and corresponding clinical information from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo) database, which were used for model validation.

Osteosarcoma methylation subtype clustering analysis

Before clustering, we first performed univariate Cox regression analysis of the methylation levels (Beta values) of all CpG sites in the DNA methylation data of 83 patients. CpG sites associated with patient prognosis were retained (P value < .05). Next, unsupervised consensus clustering of these prognosis-related CpG sites was performed based on the K-means algorithm using the ConsensusClusterPlus R package.¹⁵ The t-SNE method was performed using the Rtsne package to explore the distribution of different clusters.¹⁶ Subsequently, the relationship between clustered subtypes and overall patient survival was analyzed using patient-matched survival data using the survival R package.

Subtype immune infiltration analysis

We applied 3 independent algorithms to assess the abundance of tumor microenvironment (TME) cell infiltration. The CIBERSORT R package quantified the abundance of 22 immune cell types¹⁷; the MCPcounter algorithm was used to identify 10 immune cell lineages¹⁸; and the stromal score, immune score and tumor purity were calculated for all samples by the ESTIMATE algorithm.¹⁹ In addition, we calculated the DNA methylation of tumor-infiltrating lymphocytes (MeTIL) for all samples based on the PCA method.²⁰ Finally, we compared the differences in immune infiltration abundance, immune scores and MeTIL in different subtypes.

Copy number variation and immunotherapy analysis

Whole genome gene copy number variations (CNV) profiles and chromosome arm/gene-level variations were detected by GISTIC 2.0 (http://archive.broadinstitute.org/cancer/cga/gistic).²¹ The focal length cutoff was set to 0.5, the confidence interval was set as 0.9, and all other operational parameters were set to their default values. Meanwhile, we used OncoPrint plot to show the overview of CNVs in each OS subtype. To assess the probability of response to immunotherapy in individuals with OS, the Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) algorithm was used to predict the response to immunotherapy in each patient. We further compared the differences in immunotherapy between patients of different subtypes.

Differential expression analysis and enrichment analysis

Based on the results of the subtype survival analysis, we screened for differential genes (DEGs) between poor prognosis subtypes and good prognosis subtypes. Differential expression analysis was performed using the limma package²² and the screening thresholds were set as follows: adjusted P-value less than .05 and absolute value of Log2 fold change greater than 1. Gene set enrichment analysis (GSEA) of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) in the ranked list of DEGs were performed using the clusterProfiler R package.²³ Adjusted P value <.05 was considered statistically significant.

Construction and validation of a hypoxia-related prognostic model

Based on the results of the pathway enrichment analysis, we acquired hypoxia-related genes (HRGs) from the GSEA (https://www.gsea-msigdb.org/gsea/index.jsp) website (hallmark gene sets: h.all.v7.5.1.entrez.gmt). A total of 200 HRGs were obtained. Subsequently, univariate Cox regression analysis was performed on all HRGs to identify hypoxia-related prognostic genes associated with survival (P value < .05). Further, we used hypoxia-related prognostic genes to establish a multivariate Cox regression model, and based on the bidirectional stepwise regression method, we finally established a hypoxia-related prognostic model (HRPM). Additionally, Risk scores were calculated for each patient by using the formula of HRPM. All patients were then divided into a high-risk group and a low-risk group based on the median risk score. Survival plot, time-dependent receiver operating characteristic (ROC) curve, and patient risk heat map were plotted separately for HRPM. The GSE21257 dataset was used to validate the HRPM. Furthermore, we performed a multifactorial independent prognostic analysis of risk scores and clinical information including gender, age, and metastatic status. Finally, we compared the differences between high and low risk groups in terms of immune cell infiltration and immune score.

Results

Tumor methylation subtype clustering analysis

From a total of 385291 methylation sites, we identified 33519 potential prognostic methylation sites. As shown in Figure 1B, it shows a graph of the consensus cumulative distribution function (CDF) plot for each cluster to help the user determine the choice of the number of subtypes when the CDF reaches a maximum, at which point the consistency and clustering confidence is maximized. In addition, Figure 1C is used to represent the relative increase in cluster stability. Ultimately, based on the prognostic methylation sites, we classified the tumors into 3 subtypes (C1-3) with the help of a consensus clustering method (Figure 1A). The similarity distribution of the sample subtypes in the clusters was evaluated by visualization in the t-Distributed Stochastic Neighbor Embedding (tSNE) plot. tSNE analysis showed good separation of samples between methylation subtypes (Figure 1D). Moreover, in patients with osteosarcoma, methylation clustering subgroups are significantly associated with overall survival, according to a survival analysis. The C3 subtype patients had a much poorer overall survival rate than the other 2 subtypes (Figure 1E).

Figure 1.

Osteosarcoma clustering subtype analysis. (A) Heat map of clustering of prognostic methylation sites. (B) Consensus cumulative distribution Function (CDF) Plot. (C) Relative change in area under CDF curve for k = 2-6. (D) TSNE plot shows the distribution of the 3 clustered subtypes. (E) Each DNA methylation subtype’s survival curves.

The tumor microenvironment and immune cell infiltration

The immune cell components infiltrating the OS samples were analyzed using the CIBERSORT method and the MCPcounter algorithm, respectively (Figure 2A and C). Meanwhile, in the TARGET cohort, differences in the abundance of 3 subtypes of immune cell infiltration were estimated. The infiltration abundance of dendritic cells resting and macrophages M1 was considerably different among the 3 subtypes, according to the findings based on CIBERSORT (Figure 2B). Furthermore, with the help of the MCPcounter algorithm, we found that CD8 T cells, fibroblasts and neutrophils were significantly differentially infiltrated in subtypes (Figure 2D). The C3 subtype had the lowest stromal score and the highest tumor purity, according to the ESTIMATE algorithm, whereas the immune score did not differ substantially among the 3 subtypes (Figure 3B, C, and D). We derived MeTIL score using the methylation signature of tumor-infiltrating lymphocytes, and integrated it with immune checkpoint and the immune microenvironment to map the tumor subtype immune landscape (Figure 3A). Additionally, the C3 subtype had the lowest MeTIL score of all the subtypes (Figure 3E).

Figure 2.

Analysis of immune cell infiltration in osteosarcoma subtypes. (A) CIBERSORT shows the infiltration proportion of 22 immune cell subpopulations. (B) MCPcounter shows the infiltration proportion of 10 immune cell subpopulations. (C) Differences in immune infiltration of 22 immune cells types in 3 tumor subtypes estimated using the CIBERSORT algorithm. (D) Estimation of differences in immune infiltration of 10 immune cell types in 3 tumor subtypes using the MCPcounter algorithm. ns: no significant; *P < .05. **P < .01. ***P < .001.

Figure 3.

Comparison of immune microenvironment among osteosarcoma subtypes. (A) Immune landscape of osteosarcoma subtypes. (B) Differences in immune scores among osteosarcoma subtypes. (C) Differences in stromal scores among osteosarcoma subtypes. (D) Differences in tumor purity among osteosarcoma subtypes. (E) Differences in MeTIL scores among osteosarcoma subtypes. ns: no significant; *P < .05. **P < .01. ***P < .001.

CNV and TIDE analysis

Raw CNV data from 83 samples were analyzed with GISTIC 2.0, and the locations of significant amplifications and deletions were visualized for each chromosome (Figure 4A and B). Simultaneously, we plotted the gistic scores of overall copy number on the chromosomes (Figure 4C). In addition, we found that 14q11.2 had the highest amplification frequency, while 13q14.2 had the highest deletion frequency (Figure 4D). The TIDE algorithm was used to predict the responsiveness to immunotherapy in 3 subtypes of osteosarcoma. The results showed that C3 subtypes responded better to immunotherapy than C1 and C2 (Figure 5A). Consistent with the expected results, C3 had lower TIDE scores, dysfunction scores and exclusion scores than the other subtypes (Figure 5B, C, and D). Low TIDE scores indicate low tumor immune dysfunction, low immune escape capacity and therefore high immunotherapy response.

Figure 4.

Overall copy number variation landscape in patients with osteosarcoma. (A) Copy number amplification events in osteosarcoma gene expression levels. (B) Copy number deletion events in osteosarcoma gene expression levels. (C) The gistic scores of overall copy number on the chromosomes. (D) Percentage of copy number variants in the corresponding patient cohort (Top 20).

Figure 5.

Analysis of immunotherapy for osteosarcoma subtypes. (A) Response rate of immunotherapy among 3 subtypes of osteosarcoma. (B) Differences in TIDE scores among 3 subtypes of osteosarcoma. (C) Differences in dysfunction scores among 3 subtypes of osteosarcoma. (D) Differences in exclusion scores among 3 subtypes of osteosarcoma. ns, no significant; *P < .05. **P < .01. ***P < .001.

Differential gene screening and analysis

The above results suggest that the C3 subtype has a poorer prognosis and a better response to immunotherapy, thus we compared it to other subtypes (C1 and C2) and screened for DEGs. A total of 249 DEGs were identified, with 146 genes that were down-regulated and 103 genes that were up-regulated (Figure 6A). Following that, we used GSEA to conduct gene ontology and metabolic pathway enrichment analyses of the DEGs. The GO enrichment analysis results for biological processes showed that GO items such as translational initiation (GO:0006413), extracellular matrix organization (GO:0030198) and protein targeting to ER (GO:0045047) were enriched (Figure 6B). In the metabolic pathway enrichment analysis, we found that the hypoxia pathway (NES: 1.65) was activated while the oxidative phosphorylation pathway (NES: −1.78) was inhibited (Figure 6C).

Figure 6.

Screening and enrichment analysis of differentially expressed genes among osteosarcoma subtypes. (A) Heat map of the differential gene expression patterns. (B) GO enrichment analysis of differential gene expression. (C) GSEA enrichment results for metabolic pathways.

Hypoxia-related prognostic model construction and evaluation

A total of 29 hypoxia-related prognostic genes were screened from 200 hypoxia-related genes by univariate COX regression. Subsequently, these genes were used to construct a multivariate COX regression model, and a HRPM including 12 genes was finally obtained by a bidirectional stepwise regression method (Figure 7A). The risk score for the HRPM can be calculated using the following formula: risk score = ANXA2 * 0.203 + CASP6 * 0.366 + CAVIN1 * 2.486 + DCN * 0.671 + FAM162A * 4.871 + HMOX1 * 0.674 + MAFF * 2.098 + PDK1 * 0.282 + RBPJ * 0.431 + SDC3 * 0.586 + STC2 * 1.463 + TES * 0.566. Based on the median risk score, patients were divided into high-risk and low-risk groups. Survival risk heat map showed that higher mortality in OS patients with higher risk scores (Figure 7B). The survival analysis revealed that the mortality rate in the high-risk score group was considerably greater than in the low-risk score group (Figure 7C). The AUCs for 1-year, 2-year, and 3-year of their time-dependent ROC curves were 0.939, 0.931, and 0.909, respectively (Figure 7E). The HPRM was subsequently validated in the GSE21257 cohort, and we observed that the model performed well in the validation set (Figure 7D and F). Furthermore, a multivariate regression model was built using the patient’s clinical information and risk score, and the result indicate that risk score was an independent prognosis factor (Figure 8A).

Figure 7.

Construction and validation of a hypoxia-related prognostic model. (A) Heat map of the differential gene expression patterns. (B) GO enrichment analysis of differential gene expression. (C) Survival curves for different risk score groups in the training set. (D) Survival curves for different risk score groups in the validation set. (E) Time-dependent overall survival ROC for training set (TARGET). (F) Time-dependent overall survival ROC for validation set (GSE21257).

Figure 8.

Independent prognostic analysis of model risk scores and immune infiltration analysis. (A) Independent prognostic analysis of risk scores. (B) Infiltration abundance of 10 immune cell types in groups with high and low risk scores. (C) Comparison of immune scores in high and low risk score groups. (D) Comparison of stromal scores in high and low risk score groups. (E) Comparison of tumor purity in high and low risk score groups. ns, no significant; *P < .05. **P < .01. ***P < .001.

Differences in immune infiltration between high and low risk groups

The MCPcounter algorithm was used to calculate differences in infiltration of 10 immune cell types in high and low risk score groups. In the high-risk score group, CD8 T cells, fibroblasts, monocytic lineage, and T cells were significantly lower abundantly infiltrated than in the low-risk score group (Figure 8B). The immune score and stromal score in the low-risk group were considerably higher than in the high-risk group (Figure 8C and D), suggesting that tumor cells in the low-risk group had higher immune cell infiltration than tumor cells in the high-risk group. Correspondingly, tumor cell purity was higher in the high-risk score group than in the low-risk score group (Figure 8E).

Discussion

Previous studies have shown that OS is closely associated with aberrant genetic and epigenetic changes, leading to abnormal expression of oncogenes or methylation of antioncogenes.²⁴ Meanwhile, a large number of studies have confirmed the reliability of clustering subtype analysis of solid tumors based on methylation sites.^25-27 In this study, 33519 methylation sites were found to be statistically associated with prognosis in OS patients. Subsequently, we categorized OS into 3 clustering subtypes based on prognosis-related methylation sites, with the C3 subtype possessing the poorest prognosis of the 3 subtypes, according to the survival analysis. Meanwhile, we observed that the C3 subtype had the lowest CD8 T cell infiltration, but the highest macrophage M1 infiltration. Immune infiltration in the tumor serves a crucial function in tumor progression.²⁸ Previous studies have linked an increase in the number of CD8 T cells infiltration tumors to a favorable anti-tumor immune response, suggesting that it could be used as a prognostic indicator,^29,30 and that CD8 T lymphocyte infiltration improves the prognosis of patients with osteosarcoma.³¹ Macrophages, which are polarized M1 and M2 macrophages, are a prominent component of the tumor microenvironment.³² The formation of the tumor microenvironment encourages M1 to M2 conversion. Early tumor tissue is dominated by M1 macrophages, which limit angiogenesis and stimulate tumor immunity.³³ Contrary to expectations, we discovered that M1 macrophages of the C3 subtype exhibited the highest infiltration abundance among the 3 subtypes. In osteosarcoma, Buddingh et al. discovered that macrophages’ anticancer efficacy outweighed their probable tumor-supporting role.³⁴ Tumor-infiltrating lymphocytes (TILs) have been shown to have a profound impact on the prognosis of a variety of cancers.^35-37 Previous studies have utilized the methylation signature of TILs compared to gene expression-based immune markers to measure TILs distribution and to predict survival and tumor immune response.³⁸ Therefore, we evaluated the differences in MeTIL scores among the 3 osteosarcoma subtypes. We observed significantly lower scores for the C3 subtype than for the other 2 subtypes. In general, elevated levels of tumor-infiltrating lymphocytes are thought to be associated with a better prognosis. In particular, the synergistic effect of tumor stromal and TILs has an enhanced prognostic impact.^39,40 Based on the TIDE website, we evaluated the efficacy of immunotherapy against 3 OS subtypes, and the results suggested that the C3 subtype responded well to immunotherapy. TIDE has shown potential in published clinical trials for predicting patient response to immunotherapy, but more studies are needed to prove its utility in OS.

Next, we screened for genes that were differentially expressed between the C3 subtype and other subtypes. The pathway enrichment analysis suggested that these DEGs were substantially related with hypoxia pathway activation and oxidative phosphorylation pathway inhibition. Studies have shown that patients with lower oxygen levels are not only more tolerant to radiation than patients with higher oxygen levels, but also have a greater likelihood of recurrence.⁴¹ Hypoxia is also linked to chemotherapeutic resistance.^42,43 Anti-angiogenic medications cut off the tumor’s blood supply, causing or exacerbating a hypoxic microenvironment in the tumor, which is the primary cause of resistance to anti-angiogenic therapy after the first response, resulting in treatment failure.⁴⁴ This treatment-induced hypoxia has been proven to be the catalyst for secondary anti-vegf therapeutic resistance.^45,46 Thus, we speculate that patients with C3 subtype have a worse prognosis, which could be connected to chemotherapy resistance due to increased hypoxia pathway activation.

We obtained 200 hypoxia pathway genes from the MSigDB database, and screened 12 genes by univariate and multivariate COX regression analysis to construct HRPM and calculate the risk scores of the patients. According to the median risk score, patients were separated into high and low risk groups, and we found that patients in the high risk group had a lower overall survival rate. Furthermore, the ROC analysis revealed that the risk score model is highly reliable. We externally validated the model in the GSE21257 dataset and also obtained more satisfactory results. In a multivariate regression analysis, we included age, sex, metastasis, and risk score, and the result indicated that risk score could be an independent prognostic factor. Notably, CD8 T cell infiltration was markedly lower in individuals with high risk scores than in patients with low risk scores. Furthermore, patients in the high-risk group had lower immune score and higher tumor purity than those in the low-risk group.

Overall, in this study, we proposed a novel 12-gene signature for OS prognosis encompassing ANXA2, CASP6, CAVIN1, DCN, FAM162A, HMOX1, MAFF, PDK1, RBPJ, SDC3, STC2, and TES, which can be used to predict high-risk groups of patients with osteosarcoma. The contribution of HRPGs as biomarkers of tumor progression as well as potential therapeutic targets might be confirmed by further studies.

Conclusions

Overall, we constructed a hypoxia-associated gene signature based on multi-omics data to assess the prognosis of patients with OS. It may serve as a potential prognostic marker for the prognosis diagnosis of OS patients.

Footnotes

Acknowledgements

We thank the TCGA, GEO, and TARGET databases.

Funding:

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Declaration of Conflicting Interests:

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Declaration of Figures’ Authenticity

All figures submitted have been created by the authors who confirm that the images are original with no duplication and have not been previously published in whole or in part. All of the figures were drawn by the first author, Ye Tong.

ORCID iD

Ye Tong

References

Leichter

Sullivan

Eccle

, et al. MicroRNA expression patterns and signalling pathways in the development and progression of childhood solid tumours. Mol Cancer 2017;16:15.

Zhang

Wang

, et al. Self-stabilized hyaluronate nanogel for intracellular codelivery of doxorubicin and cisplatin to osteosarcoma. Adv Sci 2018;5: 1700821.

Zhang

Dong

Chen

, et al. Extracellular vesicle-mediated delivery of miR-101 inhibits lung metastasis in osteosarcoma. Theranostics 2020;10:411-425.

Wang

Qin

Liang

, et al. Targeting the CK1α/CBX4 axis for metastasis in osteosarcoma. Nat Commun 2020;11:1141. doi:10.1038/s41467-020-14870-4

Qiao

Yang

Xie

, et al. Engineered algae: a novel oxygen-generating system for effective treatment of hypoxic cancer. Sci Adv 2020;6:eaba5996.

Poljsak

Kovac

Dahmane

, et al. Cancer etiology: a metabolic disease originating from life’s major evolutionary transition? Oxid Med Cell Longev 2019;2019:1-16.

Han

C-F

Cai

, et al. The critical role of AMPK in driving Akt activation under stress, tumorigenesis and drug resistance. Nat Commun 2018;9:4728.

Alam

Tang

Maitituoheti

, et al. KMT2D deficiency impairs super-enhancers to confer a glycolytic vulnerability in lung cancer. Cancer Cell 2020; 37:599-617.e7.

Zhu

Wang

, et al. Hypoxia-induced TET1 facilitates trophoblast cell migration and invasion through HIF1α signaling pathway. Sci Rep 2017;7:8077.

10.

Zhao

Chen

, et al. The DNA methylation landscape of advanced prostate cancer. Nat Genet 2020;52:778-789.

11.

Deng

Yuan

Ren

, et al. A four-methylated LncRNA signature predicts survival of osteosarcoma patients based on machine learning. Genomics 2021;113: 785-794.

12.

Rosenblum

Wijetunga

Fazzari

, et al. Predictive properties of DNA methylation patterns in primary tumor samples for osteosarcoma relapse status. Epigenetics 2015;10:31-39.

13.

Chung

H-H

Kuo

C-C

Hsiao

C-W

, et al. A novel prognostic DNA methylation panel for colorectal cancer. Int J Mol Sci 2019;20:4672.

14.

Lin

Liu

, et al. Genome-wide DNA methylation profiling identifies epigenetic signatures of gastric cardiac intestinal metaplasia. J Transl Med 2020; 18:292.

15.

Wilkerson

Hayes

DN.

ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-1573.

16.

Van Der Maaten

. Accelerating t-SNE using tree-based algorithms. J Mach Learn Res 2014;15:3221-3245.

17.

Newman

Liu

Green

, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-457.

18.

Becht

Giraldo

Lacroix

, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2019;17:218.

19.

Yoshihara

Shahmoradgoli

Martínez

, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4: 2612.

20.

Jeschke

Bizet

Desmedt

, et al. DNA methylation–based immune response signature improves patient diagnosis in multiple cancers. J Clin Invest 2017127: 3090-3102.

21.

Mermel

Schumacher

Hill

, et al. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol 2011;12:R41.

22.

Ritchie

Phipson

, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43:e47-e47.

23.

Wang

L-G

Han

Q-Y.

clusterProﬁler: an R Package for comparing biological themes among gene clusters. OMICS J Integr Biol 2012;16:284-287.

24.

, et al. MicroRNA-33b suppresses migration and invasion by targeting c-Myc in osteosarcoma cells. PLoS ONE 2014;9:e115300.

25.

Wang

Zhang

Xie

, et al. Detecting the long non‑coding RNA signature related to spinal cord ependymal tumor subtype using a genome‑wide methylome analysis approach. Mol Med Rep 2019;20:1531-1540.

26.

Hill

Kelly

Kuijjer

, et al. An imprinted non-coding genomic cluster at 14q32 defines clinically relevant molecular subtypes in osteosarcoma across multiple independent datasets. J Hematol Oncol 2017;10:107.

27.

O’Sullivan

Johnson

Skinner

, et al. Epigenetic and genetic burden measures are associated with tumor characteristics in invasive breast carcinoma. Epigenetics 2015;11:344-353.

28.

Arnold

Opdenaker

Flynn

Sims-Mourtada

Wound healing and cancer stem cells: inflammation as a driver of treatment resistance in breast cancer. Cancer Growth Metastasis 2015;8:CGM.S11286.

29.

Yoon

J-Y

Gedye

Paterson

Ailles

Stem/progenitor cell marker expression in clear cell renal cell carcinoma: a potential relationship with the immune microenvironment to be explored. BMC Cancer 2020;20:272.

30.

Vidotto

Nersesian

Graham

, et al. DNA damage repair gene mutations and their association with tumor immune regulatory gene expression in muscle invasive bladder cancer subtypes. J Immunother Cancer 2019;7:148.

31.

Fritzsching

Fellenberg

Moskovszky

, et al. CD8 ⁺ /FOXP3 ⁺ -ratio in osteosarcoma microenvironment separates survivors from non-survivors: a multicenter validated retrospective study. OncoImmunology 2015;4:e990800.

32.

Marinari

Allard

Gustave

, et al. Inflammation and lymphocyte infiltration are associated with shorter survival in patients with high-grade glioma. OncoImmunology 2020;9:1779990.

33.

Han

Yang

Sheng

, et al. The mechanism of lncRNA-CRNDE in regulating tumour-associated macrophage M2 polarization and promoting tumour angiogenesis. J Cell Mol Med 2021;25:4235-4247.

34.

Buddingh

Kuijjer

Duim

RAJ

, et al. Tumor-infiltrating macrophages are associated with metastasis suppression in high-grade osteosarcoma: a rationale for treatment with macrophage activating agents. Clin Cancer Res 2011;17:2110-2119.

35.

Chen

Liu

, et al. Immune modulation of liver sinusoidal endothelial cells by melittin nanoparticles suppresses liver metastasis. Nat Commun 2019;10:574.

36.

Yang

Huang

The tumor immune microenvironment and immune-related signature predict the chemotherapy response in patients with osteosarcoma. BMC Cancer 2021;21:581.

37.

Thibaudin

Chaix

Boidot

, et al. Human ectonucleotidase-expressing CD25 ^high Th17 cells accumulate in breast cancer tumors and exert immunosuppressive functions. OncoImmunology 2016;5:e1055444.

38.

Guo

Zhu

, et al. A four-DNA methylation biomarker is a superior predictor of survival of patients with cutaneous melanoma. eLife 2019;8:e44310.

39.

Yang

Rao

Luo

, et al. Relapse-related molecular signature in early-stage lung adenocarcinomas based on base excision repair, stimulator of interferon genes pathway and tumor-infiltrating lymphocytes. Cancer Sci 2020;111:3493-3502.

40.

da Silva

Yoshida

Cardozo

, et al. Natural killer cells response to IL-2 stimulation is distinct between ascites with the presence or absence of malignant cells in ovarian cancer patients. Int J Mol Sci 2017;18:856.

41.

Shahrzad

Bertrand

Minhas

Coomber

Induction of DNA hypomethylation by tumor hypoxia. Epigenetics 2007;2:119-125.

42.

Cosse

J-P

Michiels

Tumour hypoxia affects the responsiveness of cancer cells to chemotherapy and promotes cancer progression. Anticancer Agents Med Chem 2008;8:790-797.

43.

O’Reilly

Gubbins

Sharma

, et al. The fate of chemoresistance in triple negative breast cancer (TNBC). BBA Clin 2015;3:257-275.

44.

Bergers

Hanahan

Modes of resistance to anti-angiogenic therapy. Nat Rev Cancer 2008;8:592-603.

45.

Loges

Schmidt

Carmeliet

Mechanisms of resistance to anti-angiogenic therapy and development of third-generation anti-angiogenic drug candidates. Genes Cancer 2010;1:12-25.

46.

Mavroeidis

Sheldon

Briasoulis

, et al. Metronomic vinorelbine: anti-angiogenic activity in vitro in normoxic and severe hypoxic conditions, and severe hypoxia-induced resistance to its anti-proliferative effect with reversal by Akt inhibition. Int J Oncol 2015;47:455-464.