A novel signature derived from immunoregulatory and hypoxia genes predicts prognosis in liver and five other cancers

Abstract

Background

Despite much progress in cancer research, its incidence and mortality continue to rise. A robust biomarker that would predict tumor behavior is highly desirable and could improve patient treatment and prognosis.

Methods

In a retrospective bioinformatics analysis involving patients with liver cancer (n = 839), we developed a prognostic signature consisting of 45 genes associated with tumor-infiltrating lymphocytes and cellular responses to hypoxia. From this gene set, we were able to identify a second prognostic signature comprised of 8 genes. Its performance was further validated in five other cancers: head and neck (n = 520), renal papillary cell (n = 290), lung (n = 515), pancreas (n = 178) and endometrial (n = 370).

Results

The 45-gene signature predicted overall survival in three liver cancer cohorts: hazard ratio (HR) = 1.82, P = 0.006; HR = 1.84, P = 0.008 and HR = 2.67, P = 0.003. Additionally, the reduced 8-gene signature was sufficient and effective in predicting survival in liver and five other cancers: liver (HR = 2.36, P = 0.0003; HR = 2.43, P = 0.0002 and HR = 3.45, P = 0.0007), head and neck (HR = 1.64, P = 0.004), renal papillary cell (HR = 2.31, P = 0.04), lung (HR = 1.45, P = 0.03), pancreas (HR = 1.96, P = 0.006) and endometrial (HR = 2.33, P = 0.003). Receiver operating characteristic analyses demonstrated both signatures superior performance over current tumor staging parameters. Multivariate Cox regression analyses revealed that both 45-gene and 8-gene signatures were independent of other clinicopathological features in these cancers. Combining the gene signatures with somatic mutation profiles increased their prognostic ability.

Conclusions

This study, to our knowledge, is the first to identify a gene signature uniting both tumor hypoxia and lymphocytic infiltration as a prognostic determinant in six cancer types (n = 2712). The 8-gene signature can be used for patient risk stratification by incorporating hypoxia information to aid clinical decision making.

Background

Hepatocellular carcinoma (HCC) is the sixth most common cancer and the second leading cause of cancer-related mortality worldwide [1]. Due to an initially asymptomatic disease course, this aggressive cancer once diagnosed has especially poor outcomes. The etiological risk factors for HCC varies across geographical locations [2]. This pattern mirrors the burden of viral hepatitis in Asia and Africa. Whereas in the Western world, the risk can be attributed more to alcoholic and non-alcoholic steatohepatitis. Therapy for HCC, involving either surgery or radiological ablation, is most effective when the cancer is detected early. Prompt diagnosis, however, requires regular liver imaging, usually 6-monthly ultrasound scans, which is both resource intensive and dependent upon identification of otherwise silent risk factors. Once diagnosed, curative liver transplantation is based on tumor size and restricted to patients where the primary cancer is thought unlikely to recur in the new liver. The use of size as the only prognostic criterion precludes therapy in patients with large, indolent cancers that have a low probability of recurrence and many centers are questioning whether additional data on tumor biology would add value to the diagnostic pathway.

Gene signatures obtained from tumor that are derived from common oncogenic pathways can be used to risk-stratify patients to provide individualized care. Despite genomic instability driving intertumoral heterogeneity, solid malignancies do share two common characteristics. Solid tumors are often hypoxic due to aberrant vasculature resulting in metabolic shifts towards aerobic glycolysis known as the Warburg effect [3, 4]. Tumor hypoxia is associated with metastases and aggressive phenotypes that influence clinical outcomes [5, 6]. Additionally, within tumor cells one can find lymphocytes termed tumor-infiltrating lymphocytes. These can affect patient prognosis in multiple cancers [7, 8]. A subset of lymphocytic infiltrates known as the FoxP3⁺ regulatory T cells (Tregs) function to suppress cytotoxic T cells activity to maintain tumor tolerance. Increases in Tregs contributes to unfavorable prognosis in multiple solid tumors [9]. Importantly, infiltrating Tregs within the tumor milieu can be affected by hypoxia where the latter facilitates Tregs recruitment to promote angiogenesis and maintain growth [10]. Another category of T cells known as the cytotoxic CD8⁺ T cells are required to kill cancer cells. But they can be made ineffective by Tregs or be rendered dysfunctional when exhausted [8]. Together, accumulation of Tregs and exhausted CD8⁺ T cells further exacerbate the immunosuppressive functions in solid tumors.

By taking advantage of the intricate link between hypoxia and Treg infiltration in solid malignancies, we systematically developed two prognostic gene signatures for risk stratification. Both signatures can accurately predict high-risk patients as confirmed by multi-cohort validations in six cancer types to support their validity and immediate clinical application.

Methods

Study cohorts

HCC data sets (n = 839) used in this study were GSE14520 (n = 242) from the Liver Cancer Institute, LIRI-JP (n = 226) from the International Cancer Genome Consortium (ICGC) database and TCGA-LIHC (n = 371) from the Cancer Genome Atlas (TCGA) Research Network (Additional file 1) [11]. Full clinical characteristics of HCC patients were listed in Additional file 2. Among these patients, underlying etiology differed between the cohorts with GSE14520 comprising of patients with HBV, while LIRI-JP and TCGA-LIHC consisting of patients with HBV, HCV and/or non-alcoholic steatohepatitis.

Gene expression data for 24 other cancer types (n = 20,241) generated by TCGA Research Network [11] were downloaded from Broad Institute GDAC Firehose. Illumina HiSeq rnaseqv2 Level 3 RSEM normalized gene expression profiles were converted to log₂(x + 1) scale. Expression values of genes associated with tumor-infiltrating T cells were corrected for levels of infiltration by dividing expression values of each gene with CD3D values.

Differential expression analysis

To determine differentially expressed genes between tumor and adjacent normal liver tissues in the GSE14520 cohort, the Bayes method and linear model were employed using the R package limma. P-values were adjusted using the false discovery rate controlling procedure of Benjamini–Hochberg. Genes with log₂ fold change of > 1 and adjusted P-values < 0.05 were considered significant.

Gene signature and risk scores

Expression scores for hypoxia, tumor-infiltrating T cells, 45-gene signature and 8-gene signature were calculated for each patient in all the cohorts by taking the average log₂ expression values. Nonparametric Mann–Whitney–Wilcoxon test was used to compare the distribution of hypoxia scores in tumor and adjacent non-tumor samples. Risk scores for each patient were determined by taking the sum of Cox regression coefficient for each of the signature genes multiplied with its corresponding expression value. Nonparametric Spearman’s rank correlation was employed to assess the relationship of gene scores and risk scores with tumor hypoxia (hypoxia score).

Survival analyses

We employed the Cox proportional hazards model to investigate the association between patient survival duration and one or more risk factors, e.g. 45-gene or 8-gene score, tumor stage and other clinical variables. Univariate analysis was performed to determine the influence of individual risk factors on overall survival. As multiple covariates can potentially influence patient prognosis, multivariate analyses were performed by including risk factors that were significantly associated with overall survival identified in univariate analysis (P < 0.05). Hazard ratios (HR) were determined from Cox models with HR greater than one indicating that a covariate was positively associated with event probability (increased hazard) and negatively associated with survival length. Cox regression analyses were performed using the R survival and survminer packages. Proportional hazards assumption was supported by a non-significant relationship between scaled Schoenfeld residuals and time using the R survival package. In addition, Kaplan–Meier and log-rank tests were used in univariate analyses of the gene signatures in relation to patient survival and were performed using the survival and survminer packages. Patients were median-dichotomized into low and high-risk groups based on mean expression scores of signature genes and Kaplan–Meier estimates of survival probability over time were generated. Difference between high and low-risk groups were tested using the log-rank test.

Time-dependent receiver operating characteristic (ROC) curve analysis was used to assess the predictive performance of both 45-gene and 8-gene signatures in comparison with standard tumor staging parameters. The R survcomp package was employed to compute time-dependent ROC curves. ROC curves depicted true positive rates (sensitivity) versus false positive rates (1-specificity). The area under the curve (AUC) is a measure of how well the gene signatures can predict patient survival where AUC ranges from 0.5 (for an uninformative marker) to 1 (for a perfect predictive marker).

Visualization of sample distance in the reduced 2-dimensional space

To determine whether the gene signatures (45-gene or 8-gene signature) can distinguish tumor and non-tumor samples, multidimensional scaling (MDS) analyses were performed using the R vegan package to visualize samples’ distance in the reduced 2-dimensional space. Euclidean genetic distances between each sample were investigated by MDS ordination. Permutational multivariate analysis of variance (PERMANOVA) was employed to test for differences between tumor and non-tumor samples.

Biological enrichment and protein interaction network analyses

Analysis of biological pathway enrichment on the 45-gene set was conducted using GeneCodis [12] against the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases. The Enrichr [13] tool was used to identify transcription factors from the ChEA database that are potential regulators of the 45 genes. Functional protein association network of the 45 genes was determined using the STRING [14] database.

Somatic mutation identification

Level 3 mutation datasets were downloaded from GDAC. Kaplan–Meier analysis and log-rank tests were employed to determine the association of somatic mutations in combination with the 45-gene or 8-gene signatures, on overall survival.

All graphs were generated using the ggplot2 package in R. Heatmap was generated using the R pheatmap package.

Results

Hypoxia is associated with immunosuppression through enhanced expression of Tregs and exhausted CD8⁺ T cells genes

Since both hypoxia and high density of infiltrating Tregs are linked to poor clinical outcomes, we hypothesized that these features could be employed to predict prognosis in patients with hepatocellular carcinoma (HCC). Hypoxia inducible factors (HIFs) are transcription factors (TFs) that play key roles in regulating cellular responses to hypoxia [15]. The alpha subunits (HIF-1α and -2α) heterodimerize with the beta subunit (HIF-1β) to orchestrate physiological responses in hypoxia [16]. To develop our gene signature uniting hypoxia and tumor-infiltrating T-cells, we utilized three gene sets: (1) pan-cancer hypoxia genes (52 genes) [6]; (2) HCC-infiltrating Tregs (401 genes) and (3) HCC-infiltrating exhausted CD8⁺ T cells (82 genes) [8] (Fig. 1). We utilized a HepG2 hepatoma HIF chromatin immunoprecipitation sequencing (ChIP-seq) dataset, GSE120885 [17], to determine which of these genes were bound by HIFs (Fig. 1).

Fig. 1

Schematic diagram of the study design and development of gene signatures. A liver cancer cohort (GSE14520) was used to define the first 45-gene signature. Briefly, 79 tumor-infiltrating T-cell genes were identified as HIF targets using a HIF-1α/2α ChIP-seq dataset. Of these 79 genes, 26 genes were > 1.5-fold upregulated in GSE14520. Independently, 23 hypoxia genes were identified as HIF targets. Uniting the 23 hypoxia-HIF genes and 26 T-cell-HIF genes resulted in 45 unique genes representing the first gene signature. Cox regression analyses of individual 45 genes in each of the three liver cancer cohorts (GSE14520, TCGA-LIHC and LIRI-JP) revealed a common prognostic set consisting of 8 genes that represent the second signature. This 8-gene signature is further validated in liver and five other cancers using Kaplan–Meier, Cox regression and receiver operating characteristic analyses

To determine the extent of hypoxia in cancers, we interrogated in vivo mRNA expression patterns of the 52 hypoxia genes [6] in 25 cancer types including HCC (n = 20,662) retrieved from The Cancer Genome Atlas (TCGA) [11] (Additional file 1). Hypoxia scores for each tumor and non-tumor samples were determined by obtaining the mean expression values (log₂) of the 52 hypoxia genes [6]. We observed that tumor samples were significantly more hypoxic than non-tumor samples in 20 out of 25 cancers, which included HCC (Additional file 3A). Multidimensional scaling analyses revealed that the 52 genes can distinguish tumor from non-tumor samples in these cancers, hence hypoxic transcriptional states can be used as a proxy for identifying cancerous cells (Additional file 3B).

We predict that patients with more hypoxic tumors would have higher expression of tumor-infiltrating T cell genes since hypoxia could promote the maintenance of immunological escape via the suppressive function of Tregs [10]. The HCC-infiltrating T cell gene set uniting both Tregs and exhausted CD8⁺ T cells consisted of 438 genes collectively, with 45 genes found to be in common between the two cell types (Fig. 1). Infiltrating T cell expression scores for each patient were calculated as mean expression values (log₂) of the 438 genes (Fig. 1). Patient hypoxia scores significantly correlated with tumor-infiltrating T cell scores in three independent HCC cohorts (n = 839; Additional file 2): GSE14520 (ρ = 0.564, P < 0.0001); TCGA-LIHC (ρ = 0.390, P < 0.0001) and LIRI-JP (ρ = 0.633, P < 0.0001) (Additional file 3C), suggesting that highly hypoxic tumors have transcriptional profiles associated with immunosuppressive functions.

Since HIFs are the key TFs that regulate numerous important aspects of oncogenesis [18], we hypothesize that identifying genes that are directly bound by HIFs may have more profound implications on prognosis. Using the HIF ChIP-seq dataset, we observed that 79 of the 438 tumor-infiltrating T cell genes were direct HIF targets (Fig. 1). Notably, Spearman’s correlation coefficients between the 79 genes and hypoxia score were significantly higher: GSE14520 (ρ = 0.669, P < 0.0001), TCGA-LIHC (ρ = 0.469, P < 0.0001) and LIRI-JP (ρ = 0.694, P < 0.0001) (Additional file 3D), suggesting a role for HIFs in regulating tumor lymphocytic infiltration in HCC.

Levels of T cell infiltration were corrected by dividing the expression of each of the 79 genes with CD3D expression. Of the 79 genes, 26 were upregulated more than 1.5-fold in the training cohort, GSE14520 (Additional file 3E). On the other hand, ChIP-seq analysis revealed that 23 of the 52 hypoxia genes mentioned earlier were direct HIF targets (Fig. 1). The final gene set was comprised of 45 genes that encompass 26 tumor-infiltrating lymphocyte genes and 23 hypoxia genes, with 4 genes being in common (Fig. 1; Additional file 4).

The 45-gene signature is strongly associated with poor overall survival in HCC

We next assessed the ability of the 45-gene signature (HCC45) to predict overall survival (OS) in three independent HCC cohorts (n = 839). The signature can successfully discriminate between tumor and non-tumor samples in these cohorts (Additional file 5A). Tumor samples were notably less tightly clustered implying significant intertumoral heterogeneity (Additional file 5A). We next determined HCC45 scores, calculated for each patient as the mean expression of the 45 genes. Patients were ranked according to their HCC45 scores and divided into high- or low-risk groups using the median cutoff. The OS rates were significantly reduced in high-risk patient groups with log-rank P values of 0.012, 0.0042 and 0.0043 in the GSE14520, TCGA-LIHC and LIRI-JP cohorts respectively (Fig. 2a). Disease-free survival rates were also significantly lower in high-risk patients (P = 0.026) (Additional file 5B).

Fig. 2

Patient stratification using the 45-gene signature in HCC cohorts. a Kaplan–Meier plots of overall survival in HCC patients across three cohorts stratified into low- and high-risk groups using the 45-gene signature. P-values were calculated from the log-rank test. b Kaplan–Meier plots show independence of the signature over current staging systems in HCC cohorts. Patients were sub-grouped according to TNM stages and further stratified using the 45-gene signature. The signature successfully identified high-risk patients in different TNM stages. P-values were calculated from the log-rank test. c Analysis of specificity and sensitivity of the signature in HCC cohorts with receiver operating characteristic (ROC). Plots depict comparison of ROC curves of signature and clinical tumor staging parameters. The signature demonstrates an incremental value over current staging systems. AUC: area under the curve. TNM tumor, node, metastasis staging. d Kaplan–Meier plots depicting combined relation of TP53 or CTNNB1 mutation status with the signature on overall survival in HCC

To test the independence of HCC45 over current staging systems, we performed subgroup analysis by applying the signature to patients within each tumor, node, metastasis (TNM) stages. The signature successfully identified high-risk patients, particularly among early-stage (1 and/or 2) patients (Fig. 2b). To evaluate the predictive performance of HCC45 on 5-year OS, we performed the receiver operating characteristic (ROC) analysis in comparison with tumor staging parameters. The area under the curves (AUC) for HCC45 were 0.759 (GSE14520), 0.804 (TCGA-LIHC) and 0.883 (LIRI-JP) (Fig. 2c). All HCC45 AUCs were higher than that of TNM staging (Fig. 2c). Moreover, combining HCC45 with TNM staging considerably improved its predictive ability: AUC = 0.800 (GSE14520), 0.830 (TCGA-LIHC) and 0.890 (LIRI-JP) (Fig. 2c).

We assessed two common mutations in HCC, TP53 and CTNNB1. Mutation in the tumor suppressor TP53 is associated with poor prognosis [19, 20]. Results on β‐catenin CTNNB1 mutations have been mixed; some demonstrating poor outcomes [21, 22], while others demonstrating favorable prognosis [23]. Both TP53 (P < 0.001) and CTNNB1 (P = 0.025) mutations were associated with poor prognosis (Table 1). Mortality in patients with high HCC45 scores and TP53 mutation were ~ 45% higher than in patients with low HCC45 scores harboring wild-type TP53 at 2 years (Fig. 2d). Likewise, joint relation of HCC45 and CTNNB1 predicted a ~ 35% increase in mortality at 4 years between the lowest and highest risk groups (Fig. 2d).

Table 1 Univariate and multivariate Cox proportional hazard regression analysis of risk factors associated with overall survival

We next examined the relation of HCC45 to other clinical parameters and observed that it remained highly prognostic after adjustment for tumor size, cirrhosis, fibrosis, TNM stages, Barcelona Clinic Liver Cancer stages, vascular invasion, TP53/CTNNB1 mutation status and/or alpha-fetoprotein levels in a multivariate Cox regression analysis: GSE14520 hazard ratio [HR] (HR = 1.89, P = 0.0052), TCGA-LIHC (HR = 1.61, P = 0.033) and LIRI-JP (HR = 2.20, P = 0.024) (Table 1). Both cirrhosis and HCC45 harbor complementary prognostic information contributing to elevated HR to 4.59 (95% CI 1.11–19.03, P = 0.036) (Table 1).

Gene Ontology analysis of HCC45 revealed enrichment of biological pathways linked to hypoxia, metabolism, inflammation and cancer (Additional file 6). Moreover, HCC45 was enriched for targets of several well-known cancer-related transcription factors such as MYC and RUNX1 (Additional file 6) and were functionally connected as a group (protein–protein interaction network enrichment: P < 1e−16) (Additional file 7).

The minimal prognostic 8-gene signature is sufficient and effective in predicting overall survival in HCC

Having identified that HCC45 predicts outcome in HCC, we evaluated the performance of each gene as a prognostic factor (Fig. 1). Of the 45 genes, univariate Cox regression analysis revealed that 14, 19 and 28 genes were significantly associated with unfavorable prognosis in GSE14520, TCGA-LIHC and LIRI-JP respectively (Figs. 1; 3a). Eight genes (HCC8) were found to confer prognostic information in all HCC cohorts: CA9, CCL20, CORO1C, CTSC, LDHA, NDRG1, PTP4A3 and TUBA1B (Fig. 3a; Additional file 8). Expression of HCC8 not only distinguished tumor from non-tumor sample (Additional file 9) but also increased according to tumor progression (Additional file 10). As expected, the signature successfully identified high-risk patients in full HCC cohorts, GSE14520 (P = 0.0001), TCGA-LIHC (P = 0.00016) and LIRI-JP (P = 0.0003) cohorts, and in patients stratified by tumor stage (Fig. 3b). The 5-year disease-free survival rates were also lower in high-risk patients as determined by the 8-gene signature (Additional file 11). Furthermore, as measured by Wald Chi square statistics and log-rank tests, HCC8 is superior in predicting death than HCC45 (Figs. 2a; 3b; Table 1).

Fig. 3

Minimal prognostic 8-gene signature in HCC. a Forest plots depict Cox proportional hazards analysis on 45 signature genes in three HCC cohorts. Hazard ratios (HR) were denoted as dark blue circles and light blue bars represent 95% confidence interval. Eight genes are consistently prognostic across all three cohorts, thereby constituting the minimal prognostic signature. Significant Wald test P values were indicated in red. Signature genes were highlighted in red. b The 8-gene signature successfully identified high-risk patients in different TNM stages. Kaplan–Meier plots of overall survival in HCC patients across three cohorts stratified by 8-gene signature into low and high-risk groups. Patients were stratified by the signature as a full cohort, or sub-grouped according to TNM stages. P-values were calculated from the log-rank test. Plots show independence of signature over current staging systems. c Analysis of specificity and sensitivity of the signature in HCC cohorts with ROC. Plots depict comparison of ROC curves of signature and clinical tumor staging parameters. AUC: area under the curve. TNM tumor, node, metastasis staging. d Kaplan–Meier plots depicting combined relation of TP53 or CTNNB1 mutation status with the 8-gene signature on overall survival in HCC

Multivariate models of HCC8 revealed that this reduced gene set sufficiently served as an independent prognostic risk factor: GSE14520 (HR = 1.72, P = 0.02), TCGA-LIHC (HR = 2.23, P = 0.0009) and LIRI-JP (HR = 2.58, P = 0.012) (Table 1). Predictive value for 5-year OS increased when the HCC8 signature was used in combination with TNM staging: GSE14520 AUC = 0.67 versus 0.73 and TCGA-LIHC AUC = 0.73 versus 0.76 (Fig. 3c). When TP53 or CTNNB1 mutation status and HCC8 were jointly used, patients with high HCC8 levels and mutant TP53 (P < 0.0001) or CTNNB1 (P < 0.0001) had the lowest chances of survival (Fig. 3d).

We calculated risk scores for each patient by taking the sum of Cox regression coefficient for each of the 8 genes multiplied with its corresponding expression value [24]. Risk scores significantly correlated with hypoxia scores in all 3 cohorts, indicating that patients with more hypoxic tumors have higher risk for death as predicted by HCC8 (Additional file 12), overall suggesting that it was sufficient and effective in predicting death.

Prognosis of the 8-gene signature in 5 other cancers

Hypoxic tumors recruit CD4⁺ Tregs to suppress effector T cell function and promote tumor tolerance [10]. We hypothesized that HCC8 could predict outcome in other solid malignancies in addition to HCC. HCC8 could indeed discriminate between tumor and non-tumor samples in five other cancer types (Additional file 13). Remarkably, our prognostic 8-gene signature derived from HCC was a significant adverse prognostic factor for OS in these five cancer types: head and neck (HR = 1.64, P = 0.004), renal papillary cell (HR = 2.31, P = 0.047), lung (HR = 1.45, P = 0.03), pancreas (HR = 1.96, P = 0.006) and endometrium (HR = 2.33, P = 0.003) (Fig. 4a; Additional file 14). Multivariate analysis adjusting for tumor stage revealed that the HCC8 signature retained independent prognostic relevance in these cancers: head and neck (HR = 1.66, P = 0.003), lung (HR = 1.45, P = 0.029), pancreas (HR = 1.91, P = 0.009) and endometrium (HR = 1.99, P = 0.016) (Additional file 14). Importantly, while TNM staging was not a reliable predictor of OS in pancreatic cancer on its own, our 8-gene signature successfully predicted high-risk patients when accounting for tumor stage (P = 0.009) (Additional file 14).

Fig. 4

Prognosis of the 8-gene signature in 5 other non-HCC cancers. a Kaplan–Meier plots of overall survival in patients across multiple cancers stratified into low and high-risk groups using the prognostic 8-gene signature. P-values were obtained from the log-rank test. b Analysis of specificity and sensitivity of the 8-gene signature in multiple cancers. Plots depict comparison of ROC curves of the 8-gene signature and clinical tumor staging parameters. AUC: area under the curve. TNM tumor, node, metastasis staging

The predictive value of the 8-gene signature outperformed the TNM staging system: head and neck (AUC = 0.610 versus 0.606), renal papillary cell (AUC = 0.701 versus 0.640) and pancreas (AUC = 0.694 versus 0.593) (Fig. 4b). When used in conjunction with staging information, additive effects on the predictive performance of HCC8 were observed: head and neck (AUC = 0.661), renal papillary cell (AUC = 0.768), lung (AUC = 0.671) and pancreas (AUC = 0.701) (Fig. 4b). Importantly, risk scores derived from HCC8 showed strong positive trends with tumor hypoxia, suggesting that high risk patients had more aggressive tumors: head and neck (ρ = 0.729, P < 0.0001), renal papillary cell (ρ = 0.766, P < 0.0001), lung (ρ = 0.856, P < 0.0001), pancreas (ρ = 0.844, P < 0.0001) and endometrium (ρ = 0.374, P < 0.0001) (Additional file 15).

TP53, KRAS and CDKN2A mutations are commonly observed in lung cancer [25]. CDKN2A (P = 0.03), but not TP53 (P = 0.2) or KRAS (P = 0.4), was significantly associated with poor prognosis (Additional file 14). When used in combination with our HCC8 signature, patients with high HCC8 score and mutant CDKN2A had ~ 30% higher chances for death at 2 years than patients with low HCC8 and wild-type CDKN2A (Additional file 16).

Discussion

We developed a novel gene signature that predicts overall survival in six cancers. This represents a significant step forward in the prognostic pathway since pan-cancer genetic signatures are limited at best. Capitalizing two biological phenomena of solid tumors, hypoxia and immune cells infiltration, our signature can be applied across multiple cancers, suggesting that a high degree of commonality exists in host immune response within a hypoxic tumor microenvironment. Signature genes may provide insights into biological features of individual tumors as they are implicated in oncogenic processes, i.e. pH regulation, invasion, cell proliferation, adhesion and migration, glycolysis and inflammation [26,27,28,29,30]. Combining the RNA gene signature with DNA somatic mutations simultaneously improved its prognostic capabilities, indicating that signature analysis at two macromolecular levels offers the opportunity for multiple levels of drug targeting. High-risk patients, as predicted by our signature, that concurrently have TP53 mutations suggests a role for hypoxia in promoting genomic instability and aberrant DNA damage repair.

Canonical tumor staging parameters are useful, but they do not offer additional resolution for discriminating patients with similarly-staged malignant grades. Especially for patients with stage 1 cancer, the signature allows the incorporation of hypoxia information to identify patients that are more at risk of progressing to advanced stages and developing lethal metastasis. Additionally, we were able to differentiate tumor from non-tumor samples using the signature, suggesting that it offers additional information on how transformed cells differ from normal cells. It is interesting to speculate that the signature could be used in diagnostics since non-transformed cells would have a distinct expression profile and a deviation from this profile could indicate the onset of oncogenesis.

Tumor hypoxia has wide-ranging effects causing metabolic alteration, angiogenesis, metastasis and immune suppression [31]. Significant crosstalk exists between hypoxia and antitumor immune functions where tumor hypoxia contributes to attenuated antitumor responses [10]. We observed a significant positive correlation between Treg and hypoxia gene expressions, supporting the notion that immunosuppression is higher in patients with more hypoxic tumors (Additional file 3C, D). Cancerous cells are protected through immunosuppressive functions enhanced by tumor hypoxia and together, they contribute to chemotherapy resistance [32]. Hence, to improve patient response to treatment, there is an urgent need to incorporate hypoxic tumor assessments in clinic for an effective patient-centric strategy. Through personalized therapy, patients with more hypoxic tumors would benefit from neoadjuvant treatment using hypoxia-modifying drugs to improve response to immunotherapy and chemotherapy [18].

Elevation of Tregs and exhausted CD8⁺ T cells in tumors offers another opportunity for therapeutic intervention. Exhausted CD8⁺ T cells are overrepresented among tumor-infiltrating CD8⁺ T cells [8]. These cells have reduced levels of cytotoxic markers and high expression of the PDCD1 exhaustion marker [8]. Hence, patients with high signature scores may likely have increased Treg and exhausted CD8⁺ T cell densities. These patients could benefit from therapy with the PD1 antibody as anti-PD1 treatment can rejuvenate exhausted CD8⁺ T cells [33]. When used in combination with first-line treatments, this may dramatically improve patient prognosis.

Conclusion

Although prospective validation in clinical trials is warranted, we consider our results as supporting the implications of tumor’s hypoxic and immunologic microenvironment in influencing patients’ prognosis and potentially their response to treatment. Multi-cohort validations incorporating large sample size (n = 2712) confirmed that the 8-gene signature is robust, clinically actionable and has potential to radically change how we determine prognosis and guide therapy.