3.1 Two distinct hypoxic statuses identified by consensus clustering
We performed clustering of the consensus matrix to obtain two matrices, 0 and 1, where 1 represents clustering together and 0 represents complete separation. As shown in Fig. 1A, we found that similar samples were all clustered together in the blue module. However, there was substantial noise in the blue module shown in Fig. 1B–D. Therefore, we believed that the clustering effect was optimal at k = 2, so we divided the ovarian cancer samples into two groups (cluster.A and cluster.B). The cluster.A group contained 192 patients, while the cluster.B group contained 187. Then, we illustrated the levels of differential expression of the 15 hypoxia-related genes between these two clusters using a histogram. With the exception of MRPS17, the other 14 hypoxia-related genes showed higher expression levels in the cluster.A group (Fig. 1E). PCA indicated that the 15 hypoxia-related genes could divide the two hypoxic statuses (Fig. 1F).
Differential analysis indicated that patients in the cluster.A group had increased immune cell infiltration (Fig. 1G). The Kaplan-Meier curve showed that there was no significant difference in survival outcome between the cluster.A group and the cluster.B group (Fig. 1H). Figure 2A shows that the immune checkpoint-related genes CD274, CD80, CD86, HAVCR2, IDO1, PDCD1LG2, and VTCN1 were higher expressed in the cluster.A group (Fig. 2A). Moreover, the immune activation-related genes CXCL10, CXCL9, GZMB, PRF1, and TNF had higher expression levels in the cluster.A group (Fig. 2B). Hypoxia can also induce epithelial-mesenchymal transition (EMT) of cancers [16]. In our research, we confirmed that the EMT-related genes COL4A1, PDGFRA, TGFBR2, TWIST1, VIM, and ZEB1 were higher expressed in the cluster.A group (Fig. 2C).
3.2 GO and KEGG enrichment analyses
A total of 227 hypoxia-related DEGs between the two hypoxic statuses were selected. To understand the possible enriched functions and pathways of these DEGs in OC, GO and KEGG functional enrichment analyses were performed. The GO functional enrichment analysis indicated that the genes were mainly enriched in GO:0006954 ~ inflammatory response, GO:0071347 ~ cellular response to interleukin-1, GO:0006955 ~ immune response, GO:0005125 ~ cytokine activity, and GO:0071356 ~ cellular response to tumor necrosis factor, all of which are related to inflammatory and immune responses (Fig. 3A). Moreover, the KEGG functional enrichment analysis indicated that the genes were mainly enriched in cytokine-cytokine receptor interaction, ECM-receptor interaction, NOD-like receptor signaling pathway, Jak-STAT signaling pathway, MAPK signaling pathway, and Rap1 signaling pathway (Fig. 3B).
3.3 Identification of two distinct hypoxia gene statuses
Consensus clustering was performed to based on the 227 hypoxia-related DEGs to further confirm the validity of the initial analysis. Two distinct gene clusters was obtained (gene cluster.A and gene cluster.B) (Fig. 4A–D). Figure 3E indicates that the expression levels of most of the 15 hypoxia-related genes were higher in gene cluster.A than in gene cluster.B (Fig. 4E). PCA indicated that the 227 hypoxia-related DEGs could completely distinguish the two gene clusters (Fig. 4F). Figure 4G shows that all of the other 27 of 28 immune cell types were more highly infiltrated in the gene cluster.A group, with the exception of CD56bright natural killer cells. Moreover, the Kaplan-Meier curve indicated that the survival outcome did not differ significantly between the two gene clusters (Fig. 4H). The expression levels of immune checkpoint-related genes CD274, CD80, CD86, HAVCR2, IDO1, PDCD1, PDCD1LG2, TNFRSF8, and VTCN1 were higher in the gene cluster.A group (Fig. 5A). With the exception of TBX2, all immune activation-related genes CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, and TNF were overexpressed in the cluster.A group (Fig. 5B). The EMT-related genes COL4A1, PDGFRA, TGFBR2, TWIST1, VIM, and ZEB1 were more highly expressed in the gene cluster.A group (Fig. 5C).
3.4 PCA algorithm to quantify the hypoxic status
To quantify the hypoxic status, each patient was assigned a hypoxia score using a PCA algorithm. Most of the 15 hypoxia-related genes were higher expressed in the high-score group (Fig. 6A). PCA indicated that the 227 hypoxia-related DEGs could divide the two score groups (Fig. 6B). We compared the hypoxia score between the clusters identified using either the 15 hypoxia-related genes or the 227 hypoxia-related DEGs. The hypoxia score in cluster.A or gene cluster.A was higher than that in cluster.B or gene cluster.B (Fig. 6C and D). However, the hypoxia score showed no statistically significant difference among the different stages and grades (Fig. 6E and F). The relationship among hypoxic status, hypoxia gene status, and hypoxia score was visualized in a Sankey diagram (Fig. 6G). 28 immune cell types were highly infiltrated in the high-score group (Fig. 6H). The K-M curve indicated that the survival outcome have no differ, which was in line with our expectations (Fig. 6I). Roughly similar to previous results, the expression levels of immune checkpoint-related genes (Fig. 7A), immune activation-related genes (Fig. 7B), and EMT-related genes (Fig. 7C) were higher in the high-score group. Finally, considering that heterologous proteins produced by gene mutations can lead to an increase of immune cell infiltration, we analyzed the different gene mutations between the high- and low-score groups with different immune infiltrate statuses using the “maftools” package, and visualized by forest plots (Fig. 8A) and co-oncoplots (Fig. 8B). The results indicated that the mutation frequencies of BRCA1, PRUNE2, APOB, RELN, MDN1, ZSWIM8, MGAT4C, FAT2, and VWF were higher in the high-score group, while those of HECW1 and RB1 were on the contrary.