The gene expression data and clinical information of Colorectal cancer patients (n=288) downloaded from TCGA database (https://portal.gdc.cancer.gov/) were used as discovery cohort. The CPTAC-2 prospective data set including gene expression and clinical data of Colorectal cancer patients (n=110) obtained from cBioPortal (https://www.cbioportal.org) were used as validation cohort. Gene expression of Colorectal tissues obtained from GTEx database (https://gtexportal.org/home/)was used as normal control in the downstream analysis (n=253).

We first obtained pathways and biological processes from Molecular Signatures Database (MSigDB) C2 curated gene sets on gene set enrichment analysis (GSEA) website. Then, a total of 3524 genes in these gene sets that associated with oxygen metabolism were identified as oxygen metabolism related genes in our study. Differential expressed oxygen metabolic genes (DEG) between tumor and normal samples were analyzed in discovery and validation cohort using ‘limma’ R package, respectively. Genes with FDR< 0.05 and |log₂(FoldChange)| > 1 were extracted as differentially expressed genes. The “Pheatmap” and “ggplots” package was used to plot heatmaps and volcano maps for DEGs. Venn plots of up- and down-regulated DEGs between discovery and validation cohort were achieved using a Venn online tool (https://bioinformatics.psb.ugent.be/webtools/Venn/)

DEGs significantly associated with overall survival (OS) in the entire discovery cohort were identified using univariate Cox proportional hazards regression analyses. A P-value ≤.05 was considered statistically significant. Then, we performed the least absolute shrinkage and selection operator (LASSO) penalty Cox regression analysis to eliminate genes that might overfit the model (Combined-24). Finally, we calculated risk score (RS) for each patient by a linear combination of Cox coefficient and expression of optimal prognostic DEGs identified by multivariate Cox analysis. The risk score calculation formula was as following:

Ei and Ci represented ith gene expression and corresponding coefficient value. N is the number of optimal prognostic DEGs. Patients with RS values greater than the median were defined as high-risk groups, otherwise as low-risk groups. Kaplan-Meier analysis was conducted using the “survival” and “survminer” R package. Receiver operating characteristic curve (ROC) and the “area under the ROC” (AUC) analysis were used to evaluate the performance of the prognostic model.

We used validation cohort (CPTAC) to verify the prognostic risk model. RS of each patient in validation cohort was calculated using formula mentioned above based prognostic DEGs and coefficient identified in discovery cohort. Survival and ROC analysis were used to validate the performance of prognostic risk model.

To assess the independent prognostic value of oxygen metabolic gene-based risk models in colorectal cancer, we performed both univariate and multivariate analyses of prognostic factors using Cox proportional hazards regression. Age, gender, pathological stage and TNM stage were treated as covariates. Factors with p value< 0.05 in both univariate and multivariate Cox analysis were defined as independent prognostic indicators.

We constructed the PPI network of prognostic oxygen-metabolic genes using the PathwayCommon (https://www.pathwaycommons.org/) PPI database. Analysis of functional interactions between proteins was performed in order to elucidate the potential roles of prognostic genes in colorectal cancer tumorigenic process. The PPI networks were visualized using the Cytoscape software.

The hub genes were identified using DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress and ClusteringCoefficient algrithms with Cytoscape’s plug-in cytoHubba in the PPI network. Then, we obtained all possible transcription factors (TFs) of hub gene from ChiIP seq experimental data of human samples in the ENCODE project, and identified upstream TFs that play a role in colon cancer by calculating the expression correlation between these TFs and hub genes in our tumor samples. The correlationship of gene expression between hub genes and their TFs was conducted using the spearman method. The most relevant TF-hub gene relationship was shown by scatter plots.

The workflow of this study is shown in Figure 1 . Tumor samples (n=288) obtained from the TCGA database were regarded as the discovery cohort, while tumor samples (n=102) obtained from the CPTAC project were regarded as validation cohort. We compared expression levels of 3524 oxygen metabolic genes between tumor and normal samples in discovery set and validation set, respectively. The distributions of all genes including identified DEGs according to the two dimensions of -log₁₀(FDR) and log₂(FoldChange) were displayed by volcano maps ( Figures 2A, B ). It was found that there were 262 up-regulated and 188 down-regulated genes in the discovery set, 175 up-regulated and 27 down-regulated genes in the validation set ( Figure 2C ). To obtained the more reliable prognostic gene signature, we established the prognostic model with 72 DEGs up-regulated in both the discovery set and the validation set ( Figures 1 , 2D ). We didn't obtaied reliable down-regulated DEGs which were identified in the discovery and validation cohorts Figure 2E )

Univariate Cox regression analysis revealed that 95 OM DEGs were significantly (P<.05) associated with OS in the discovery cohort. Among them, 73 DEGs were associated with good OS, while 22 DEGs were associated with bad OS.

Based on the discovery cohort, we obtained eight candidate prognostic OM genes using Lasso Cox regression analysis. Then, we acquired five optimal genes, including FLRT2 (Fibronectin Leucine Rich Transmembrane Protein 2), ATP6V0E2 (ATPase H+ Transporting V0 Subunit E2), ODC1 (Ornithine Decarboxylase 1), SEL1L3 (SEL1L Family Member 3) and VDR (Vitamin D Receptor). Four of these genes were high hazard genes, and one gene (SEL1L3) was low hazard gene, and all these genes were up-regulated DE genes ( Table 2 ). The risk score of each tumor sample was calculated as follows: risk score = (1.0232× FLRT2 _exp) + (1.0046 ×  ATP6V0E2 _exp) + (0.9806 ×  SEL1L3 _exp) + (1.0015 × ODC1 _exp) +  (1.0493 × VDR _exp).

Based on the optimized risk score threshold, all colorectal cancer patients of discovery cohort were divided into a high-risk group (n = 30) and a low-risk group (n = 252). The K-M survival analysis shown those OS times of high-risk patients were significantly longer than that of low-risk patients (p< 0.001) ( Figure 3A ). The median survival time of patients in the high-risk group was shorter than 5 years, while that of patients in the low-risk group was longer than 10 years. From the perspective of survival rate, the 1-year, 3-year and 5-year survival rates of the high-risk group were only 72%, 63% and 35% respectively, while the corresponding survival rates of the low-risk group reached 91%, 85% and 71% respectively. In addition, the AUC values of the five-OM gene prediction model were 0.753, 0.674, and 0.714 when predicting one -, three -, and five-year OS, respectively ( Figure 3B ). To find out whether all 5 prognosis OM genes are associated with advanced stages and therefore are associated with worse prognosis, we analyzed the OS time of patients with high- and low-RS from tumor stage I, II, III and IV. Results indicated that the OS time of patients with high RS was significantly shorter than that of patients with low RS in stage II and IV, which proved the prognostic effectiveness of our 5-OM gene signature in these two stages. But in stage I and III, there was no significant difference in the OS time of patients with high- and low-RS, suggesting the prognostic limitations of the model in these two stages ( Figure S1 ).

We validated the performance of the model using the validation cohort. Patients in validation cohort were divided into high- and low- risk groups based on RS threshold determined in discovery cohort. Results indicated that 12 patients and 90 patients were categorized as high- and low- risk groups, respectively. K-M survival curves were significant different between the two risk groups (p < 0.001) ( Figure 3C ) and the AUC values at 1- and 3-year were 0.974 and 0.958 in the validation cohort, respectively ( Figure 3D ). At the same time, the RSs of patients in the high-risk group were higher than those in the low-risk group, which proved that the model had a good performance in the prognosis evaluation and monitoring of colorectal cancer.

To assess whether the RS determined by five oxygen-metabolic prognostic model is an independent prognostic indicator for patients, we carried out a univariate Cox analysis to assess the impact of risk scores and clinicopathological parameters on prognosis, such as age, gender, histological type, longest dimension, pathological stage and so on. We found that the longest dimension, lymphatic invasion, pathological stage and risk score were associated with poor outcomes of prognosis in patients ( Table 3 ). Therefore, these characteristics were included in a multivariate Cox regression analysis, which indicated that age, pathological stage and the risk score estimated based the prognostic model was an independent prognostic factor for colorectal cancer ( Table 3 ). This result indicates that there is significant potential for these oxygen-metabolic genes to predict the prognosis outcome of patients with the colorectal cancer.

We investigated the regulatory relationships between TFs and prognostic genes. Firstly, we obtained upstream TFs of each gene from the ChIP-Seq experiment in the ENCODE project (https://www.encodeproject.org). Then, we analyzed expression correlation of between TFs and genes in TCGA tumor samples to validate the TF regulation in vivo. We found that 11 cancer-related TFs including CTBP2, E2F1, EP300, ETS1, FOS, JUN, MYC, RELA, STAT1, STAT3 and TCF7L2 were significantly correlated with our prognostic genes. Among them,STAT3 and MYC were significantly correlated with all prognostic genes, in which positively correlated with FLRT2, ODC1, SEL1L3 and VDR, and negatively correlated with ATP6V0E2 ( Figure 4 ).These results indicates that the prognostic genes we screened are important downstream molecules of classic cancer driver genes like STAT3 and MYC.

In order to study the possible function and mechanism of prognostic oxygen metabolic genes in colorectal cancer, we screened KEGG cancer pathway genes that can interact with prognostic genes, and then used the metascape to analyze the gene ontology (GO) and pathway enrichment of interacting genes It was found that the interaction genes were significantly enriched in GO terms including “response to inorganic substance”, “response to xenobiotic stimulus”, “response to oxidative stress”, and the enriched pathways were “transport of small molecules”, “ion channel transport”, “mineral absorption” and so on ( Figure 5A ). To further capture the relationships between these enriched terms, we constructed a network diagram using Metascape analysis. Spots represented GO terms or pathways. Larger and connected points represented the presence of more similar genes between the terms or pathways. The “Transport of small molecules” pathway contained many genes participating in Ion channel transport, while the “response to oxidative stress” gene set contained many genes participating in cell stress and inflammatory response terms or pathways ( Figure 5B ). In addition, the PPI network showed a relationship between different genes and proteins in two sub-modules ( Figure 5B ). The “Fluid shear stress and atherosclerosis” sub-module seeded by the MYC included IL1B, TP53, PIP, MYC, IL2, NFKBIA. AGT, CXCL8, BLM and MMP2, which can identify the structural components of the extracellular matrix to provide tensile strength; the “extracellular matrix organization” sub-module included SPP1, IGFBP4, GAS6, MXRA8, and SPARCL1, which play a central role in vascular biology; the “Signaling by Interleukins” sub-module seeded by the F2 included KNG1, KRT6B, PARP1, TNFRSF1A. KRT2, ZBTB16. ABCB1, BCL2, KRT1. PML, C3, TXN, KRT6C, F2, PTK2, TNF, which could enable HIF-mediated inflammatory response during cancer development ( Figure 5C ).