Dysregulated Forkhead Box (FOX) Genes Association with Survival Prognosis, Anti-tumor Immunity, and Key Targeting Drugs in Colon Adenocarcinoma

Background: Several studies have revealed that the aberrant expressions of forkhead box (FOX) genes are associated with carcinogenesis. However, the crucial biological functions of the FOX gene in colon adenocarcinoma (COAD) remain unknown. Methods: The TCGA-COAD dataset (n=328) was utilized for determining the deregulated FOX genes and their association with functional enrichment, protein-protein interaction (PPI), survival prognosis, anti-tumor immunity, cancer-associated pathways, and biological processes in COAD. In addition, we used GSE166427 (GPL13667) as a validation cohort (n=196). Molecular docking studies were applied to perform the drug interactions. Results: The FOX genes are deregulated in the COAD (Log2 FC>0.50, P<0.05), and the PPI network of FOX members is substantially related to the enrichment of cancerous signaling, immune responses, and cellular development (FDR<0.05). A worse prognosis for overall survival in COAD individuals is connected with the subgroup of FOX transcripts (P≤0.05). FOXD4, FOXH1, and FOXS1 were identified as predictive variables in the univariate and multivariate Cox regression models (P≤0.05). FOXH1 and FOXS1 are substantially linked to the deregulated immunity in COAD (R>0.20, P<0.01). Furthermore, FOXS1 expression regulates cancer-associated pathways and biological processes (P<0.05). Moreover, FOXD4, FOXH1, and FOXS1 are genetically altered and showed diagnostic efficacy in COAD. We revealed that FOXD4, FOXH1, and FOXS1 are consistently deregulated in GSE166427 (P<0.05). Finally, molecular docking revealed that FOXH1 interacted with various drugs, including belinostat, entinostat, and panobinostat. Conclusion: The FOX genes have a strong correlation with the poor prognosis for survival, tumor immunity, cancer-associated pathways, and biochemical processes that cause the pathogenesis of COAD.


Introduction
In the human genome, there are more than 50 forkhead box (FOX) genes or proteins, divided into 19 subclasses (FOXA to FOXS) based upon sequence homology inside and outside the forkhead domain. 1,2Various FOX genes regulate crucial physiological processes, including embryogenesis, cellular homeostasis, and the immune system. 3,4The functioning of the lung, kidney, immune function, and nerves is also regulated by the FOX genes. 3,4The FOX genes promote and collaborate with certain other active gene transcription and epigenetic controllers, acting as transcriptional activators, transcriptional repressors, and pioneer factors. 2 It has been reported that the FOX genes and their translated protein are involved in the development, metabolism, cancer, and aging. 2,4Cancer genesis, growth, maintenance, advancement, and medication resistance are all closely correlated with dysfunctional FOX genes. 2 Also, FOX genes are associated with substantial cancerous biological processes, including senescence, apoptosis, DNA damage repair, and survival of cancer patients. 2,5Various FOX genes are related to tumor-promoting inflammation, immune destruction, genome instability and mutation, angiogenesis, proliferative signaling, tumor invasion, metastasis, and resisting cell death. 3alignancies in several tissues, such as the lung, breast, and thyroid, have been linked to the amplification of FOX genes. 6FOX genes act downstream of several oncogenic signaling pathways, including PI3K-AKT, ERK, Wnt, β-catenin, EGFR, and NF-κB-IKKβ, associated with CRC cancerogenesis. 7The clinical results of cancer patients are strongly correlated with deregulated FOX genes.For example, Shan et al discovered that nasopharyngeal cancer patients have a poor prognosis when their FOXJ2 (OMIM number: 619162) activity is increased. 8According to Song et al, the crucial predictive factor for overall survival (OS) in those with various malignancies is FOXM1 (OMIM number: 602341). 9The FOXD1 (OMIM number: 601091) protein was significantly elevated in the primary HNSCC cohort.Significant associations existed between its abundance, cervical node metastases, and poor overall and disease-free survival (DFS). 10In colorectal cancer, the expression of FOXM1 correlates with poor prognosis. 11In addition to the clinical outcomes, FOX genes are critically associated with immunoregulatory functions.The FOX genes are associated with immunoregulation, including regulating CD4 + T cell tolerance, thymic development, macrophage differentiation, natural killer cell effector function, and T cell activation. 4The FOXOs regulate the immune system by modulating tumor and stromal cells. 12Martin et al showed that the FOXP3 (OMIM number: 300292) is an immunosuppressive marker in human malignancies. 13These studies provided crucial information about the oncogenic roles of FOX genes associated with cancer onset, development, metastasis, tumor immunity, and signaling pathways in human malignancies.
Globally, colon cancer is the second most common cause of mortality (9.4%), and the third most commonly diagnosed disease (10.0%). 14Here, we present a comprehensive bioinformatics study to evaluate the changes in FOX gene expression levels.Also, we identified the functional enrichment of deregulated FOX genes in colon adenocarcinoma (COAD).Additionally, we investigated the involvement of FOX genes in proteinprotein interactions (PPIs) and their association with the enrichment of pathways.We further examined the connection between altered FOX genes and the survival rates of COAD individuals.Univariate, multivariable, and nomogram analysis revealed COAD's independent prognostic FOX genes.Then, we identified that the prognostic FOX genes were related to the immune content, stromal content, purity of tumor, immune signatures, immune ratio, cancer-associated signaling pathways, and cancer hallmark biological processes.Finally, we observed that the prognostic FOX genes were genetically altered and effective in diagnosing COAD individuals.

Datasets
We transformed the data into a log2-base conversion after downloading the TCGA COAD dataset from https:// portal.gdc.cancer.gov.We utilized relevant information from the TCGA COAD cohort to compare the survivorship of two patient groups in the database (https://portal.gdc.cancer.gov)(n = 328, 287 tumor samples, and 41 normal samples).In addition, we investigated the transcript levels of FOX genes in multiple types of cancer using the Oncomine repository (https://www.oncomine.org/resource/login.html).To determine the genetic changes in three prognostic FOX genes, we utilized COAD (TCGA, Firehose Legacy) individuals with mutation and CNA information (n = 220) in cBioPortal (http://www.cbioportal.org/).Finally, we used GSE166427 (GPL13667, 98 tumor samples, and 98 normal samples) as a validation cohort (n = 196). 15

Investigations of FOX Gene Variations in COAD Compared to the Control
The Limma linear model, which has a range of test conditions and predictors, is specifically made for assessing complicated studies.When comparing COAD samples (n = 287) to control samples (n = 41), we used the R package "limma" to find the relevant DEGs. 16Utilizing the TCGA COAD cohort, we identified the significant FOX genes in the COAD individuals relative to the control sample.

Gene-Set Enrichment Analysis
A computational technique called gene-set enrichment analysis (GSEA) is capable of helping assess if a predefined gene collection exhibits statistically meaningful, concordant changes between two biological contexts.Using the GSEA tool, we investigated the GSEA of the candidate genes. 17We entered the gene set into the GSEA tool to find the significant GO and pathways.The Reactome pathways 18 significantly associated with the FOX genes were identified.

Protein-Protein Interaction Network Construction
Using the NetworkAnalyst application, we developed a PPI network of the DEGs. 19In the NetworkAnalyst software, we used the STRING tool 20 by selecting interactome with a medium interactome score (400) and confidence score cutoff of 900, and required experimental evidence was also determined.The node explorer module of NetworkAnalyst was used to find the ranking genes (those with a high degree of connectivity to specific other nodes) in the PPI network. 19To find KEGG pathways connected to generated PPI, we used the function explorer module of the NetworkAnalyst tool.Using the software Cytoscape 3.8.2,we illustrated the PPI networks. 21

Survival Analysis of FOX Genes in COAD
In the two groups of individuals with colon cancer, we evaluated the OS and the DFS.According to our assessment of the clinical information from the TCGA COAD cohort, the follow-up period was 147.9 months.The survival disparities between the high-expression group (HEG) and low-expression group (LEG) of FOX genes in COAD individuals were identified using the Kaplan-Meier method (HEG > median > LEG).Utilizing the R package "survival", the survival implications of significant FOX genes in the TCGA COAD database were investigated. 22We used the "coxph" function in the R package "survival" for the univariate and multivariable Cox regression assessments of variables.The nomogram was generated using the R package "rms". 22

ESTIMATE Algorithm for Quantifying Immune Score, Stromal Score, and Tumor Purity
The R-based ESTIMATE program predicts the tumor purity, immune score, and stromal score.This technique utilizes the gene expression patterns of 141 stromal genes and 141 immune genes. 23Applying relevant gene expression matrix information, the quantity of stromal cells and immune cells infiltrating into malignant cells was determined. 23Next, we contrasted the immune, stromal, and tumor purity values between the HEG and LEG of crucial FOX genes (Wilcoxon sum rank test, P < 0.05).

Single Sample Gene Set Enrichment Analysis
Instead of using a sample population, the single sample gene set enrichment analysis (ssGSEA) approach operates on a single-sample basis.For every pairing of samples and genes set in the cancer tissues, we used the extended modules of GSEA and ssGSEA to determine the enrichment levels of immune cells, cancer-associated pathways, and hallmark biological processes. 246][27][28] We identified the ssGSEA score of various immune signatures (such as CAFs, HLA genes, immune cell infiltrate genes, immune checkpoint genes, pDCs, etc), cancerous KEGG pathways (such as cell cycle, JAK-STAT, MAPK, mTOR, WNT, etc), 29 and three key hallmark biological processes (epithelial-mesenchymal transition [EMT], angiogenesis, and hypoxia).The relationships between the enrichment levels (ssGSEA scores) of immunological signatures, pathway activity, and biological processes and the levels of FOX gene expression were investigated using Spearman's correlation test.The marker gene sets with immunological signatures, pathway activity, and biological processes are listed in Table S1 (see Supplementary file 1).

Analysis of Genetic Changes in Prognostic FOX Genes
Utilizing cBioPortal, an open-access platform for assessing genetic modifications in multimodal cancer studies, we were able to determine the genetic changes in independent prognostic FOX genes.In our work, we identified the genetic of prognostic FOX genes using the COAD (TCGA, Firehose Legacy) individuals with mutation and CNA data (n = 220) in cBioPortal (http:// www.cbioportal.org/).

Diagnostic Effectiveness of Prognostic FOX Genes in COAD
To properly evaluate ROC curves and partial areas under the curve, the "pROC" R package provides an assortment of statistical models.The area under the ROC curve (AUC) was determined using the "pROC" R package to assess the capacity to discriminate between COAD and healthy samples.The graph was shown to determine the diagnostic values of prognostic FOX genes. 30The AUC values represented the diagnostic value and differences between tumor and healthy samples for each gene, and the AUC > 0.5 for a single gene was defined as the diagnostic efficacy. 31

Validation of Key Gene Expression in an Independent Dataset
As a validation cohort, we used GSE166427 (GPL13667) to confirm the expression levels of critical genes (n = 196, 98 tumor samples, and 98 normal samples). 15This database aimed to evaluate the expression profiling of COAD and normal adjacent colon cells.The platform GPL13667 is based on Affymetrix Human Genome U219 Array.We downloaded the series matrix data and utilized the t-test to identify the DEGs between colon tumors and adjacent normal colon cells.We selected the genes with the highest value of fold change for multiple probes of a single gene.

Exploration of the Drug Compound's Interaction with Key Genes and their Molecular Docking
We employed the NetworkAnalyst 19 software for extracting the chemical-protein interaction.We inputted the gene symbol of key genes into the NetworkAnalyst 19 software and identified the networks or sub-networks between the genes and chemical compounds.We utilized the Cytoscape tool to show the retrieved drug-gene interaction.Then, we used these interacting compounds for molecular docking analysis.We downloaded the protein product of FOXH1 (OMIM number: 603621) (protein database ID: 5XOC) and all other chemical compounds interacting with this gene.We prepared the protein using Discovery studio (https://3ds.com/products-services/biovia/products).First, we eliminated ligands and water molecules from the receptor proteins.Next, we prepared the ligand using PyRx (https://pyrx.sourceforge.io/).Finally, using PyRx, we conducted a molecular docking investigation (https:// pyrx.sourceforge.io/).

Statistical Analysis
Log2FC > 0.50 (absolute value) and P < 0.05 were the cutoffs we used to identify the DEGs.For choosing the significantly enriched GOs and pathways, the FDR < 0.05 was considered.When comparing the two patient groups' survival rates using Cox regression, a P ≤ 0.05 was considered significant.The Wilcoxon sum rank test was performed to contrast the two patient groups (P < 0.05).To determine the significance levels between the two factors, we used either Pearson's or Spearman's correlation test.The levels of FOX genes and the ssGSEA scores of immunological signatures, pathways, and biological processes were analyzed using Spearman's correlation test because the data were not normally distributed (P < 0.05).We used Pearson's correlation test to examine the relationships between the expression levels of FOX genes and the expression levels of other immune-marker genes and FOX genes since the data were normally distributed (P < 0.05).During the validation of key genes in the GSE166427, the P < 0.05 cutoff was established as the significant value.We used the R package "ggplot2" to present the information from our investigation.

Identifying the Differentially Expressed FOX Genes in the COAD
We analyzed the COAD's significantly differentially expressed FOX genes relative to the normal samples (Log 2 FC > 0.50, Adjusted P < 0.0001).We found that 17 FOX genes have increased expression levels in COAD (Table 1).In contrast, 15 FOX genes exhibited a downregulation trend in the COAD (Table 1).Figure 1 shows the heatmap of the genes' transcriptional value with differential expression.After passing the threshold, the expression of additional FOX genes did not change in the TCGA COAD samples.

Association of FOX Genes with Pathway Deregulation and Functional Enrichment
We utilized the GSEA tool to identify the enriched gene ontology (GO) and deregulated pathways significantly linked to the dysfunctional FOX genes (Figure 2).We identified biological processes (such as regionalization, negative regulation of DNA binding transcription factor activity, and positive regulation of developmental process) that strongly correlate to increased FOX genes (Figure 2A).The cellular components, including chromatin and chromosome, are enriched considerably with upregulated FOX genes (Figure 2A).Moreover, we found the six molecular functions (such as DNA The character "OMIM#" is the OMIM number obtained from the website https://www.omim.org/.* indicates that the OMIM number was not found on the website (https://www.omim.org/).binding bending, DNA binding transcription repressor activity, and transcription regulator activity) that strongly correlate to increased FOX genes (Figure 2A).In addition, the downregulated FOX genes enrich GO terms and signaling pathways.We identified 63 biological processes (such as positive regulation of the biosynthetic process, animal organ morphogenesis, and tube development) that are significantly associated with downregulated FOX genes (Figure 2B and Table S2).The chromatin, chromosome, and transcription regulator complex are significantly enriched cellular components related to the decreased FOX genes (Figure 2B).Moreover, we found seven enriched considerably molecular functions related to the decreased FOX genes (Figure 2B).Furthermore, we revealed the deregulated Reactome pathways associated with differentially expressed FOX genes.We found that the Reactome FOXO mediated transcription of cell cycle genes pathway is significantly enriched with upregulated FOX genes.In addition, we found seven pathways (such as AKT phosphorylates targets in the nucleus, regulation of FOXO transcriptional activity by acetylation, and regulation of localization of FOXO transcription factors) that are significantly associated with downregulated FOX genes in COAD (Figure 2C).It suggests that COAD transcriptional activity and other carcinogenic signaling mechanisms are correlated to uncontrolled FOX genes.

FOX Genes Associated with PPI Network and Correlated with Each Other
We inputted deregulated FOX genes (Table 1) into the NetworkAnalyst software to construct the PPI network.In the NetworkAnalyst software, we used the STRING tool 20 with a medium interactome score (400), and confidence score cutoff of 900, and the required experimental evidence was also selected.In the PPI, we found that the FOX genes are associated with five sub-networks.In sub-network 1, the eight FOX genes, including FOXO1, FOXO3, FOXO4, FOXM1, FOXP3, FOXG1, FOXA1, and FOXH1 are involved in the PPI network (Figure 3A).All nodes in subnetwork 1 are listed in Table S3.Interestingly, the function module of the NetworkAnalyst software identifies that sub-network 1 is substantially related to the enriched 95 KEGG pathways (Table S4).The top 30 pathways are illustrated in Figure 3B.These significant (FDR < 0.05) pathways mainly involved with cancer (such as pathways in cancer, endometrial cancer, colorectal cancer, small cell lung cancer, and central carbon metabolism in cancer), immunity (such as Th17 cell differentiation, B cell receptor signaling pathway, and T cell receptor signaling pathway), and cellular signaling and development (such as TGF-beta signaling pathway, cell cycle, and cellular senescence).It suggests that sub-network 1 is related to the regulation of tumor immunity and tumor biology in COAD.In sub-network 2, the FOXD3 interacts with four genes (Figure 3C).In other sub-networks, FOXJ1, FOXF2, and FOXK1 are involved with the PPI network (Figure 3D-F).After determining the PPI of the FOX genes, we hypothesized that the FOX genes are interconnected (Pearson's correlation test, P > 0.05).Interestingly, we found that some FOX genes correlate with each other (Figure 3G).For example, there is a substantial positive correlation between FOXP3 and FOXF1, FOXF2, FOXD3, FOXN3, FOXP2, FOXI2, FOXS1, FOXC1, and FOXO1 levels (Figure 3G).Similarly, the expression level of FOXS1 has a significant positive correlation with the expression of FOXO1, FOXI2, FOXO4, FOXP3, FOXF1, FOXF2, and FOXD3 (Figure 3G).In contrast, FOXD2 negatively correlates with the expression level of FOXC1, FOXO1, FOXD1, FOXD4, and FOXD4L1 (Figure 3G).It indicates that the FOX genes may be regulating the expression of each other and contribute to the carcinogenesis of COAD.

Dysfunctional FOX Genes Related to Poor Survival Prognosis in the COAD
We investigated the survival prognosis of dysfunctional FOX genes (shown in Table 1) using the clinical information of the TCGA COAD cohort.We discovered that patients with COAD had a worse OS prognosis when the expression level of FOXC1, FOXD1, FOXD4, FOXD4L2, FOXH1, FOXQ1, and FOXS1 is higher than the median (HEG) (Figure 4A).In contrast, the higher expression group (HEG > median) of FOXN1 has a favorable survival prognosis in COAD patients (Figure 4A).In addition, the higher expression group (HEG > median) of FOXC1 and FOXD4L2 are poorly linked with DFS time in COAD patients (Figure 4B).
Moreover, we applied univariate Cox regression analyses of the expression of seven prognostic genes (FOXC1, FOXD1, FOXD4, FOXD4L2, FOXH1, FOXQ1, and FOXS1), age, sex, weight, stage, and genome fraction altered.The univariate Cox regression analyses identified five genes (FOXC1, FOXD1, FOXD4, FOXH1, and FOXS1), stage, and fraction genome altered as significant prognostic factors (Figure 5A).Furthermore, we investigated the multivariable Cox model with the levels of FOXC1, FOXD1, FOXD4, FOXH1, and FOXS1 genes, with COAD stage (stage I and stage II versus stage III and stage IV), and fraction genome altered as the predictor variables.Remarkably, we discovered that stage, fraction genome alterations, and the levels of three genes (FOXD4, FOXH1, and FOXS1) were the independent prognostic variables (Figure 5B).Finally, the nomogram was constructed to determine the possibility of these three prognostic FOX genes and stages of the tumor with fraction genome alterations influencing the prognostic outcome (Figure 5C).

The Dysfunctional Prognostic FOX Genes Regulating the Tumor Immunity in the COAD
In human CRC, the level of tumor-infiltrating lymphocytes (TILs) is a remarkable indicator of eventual survival. 32We investigated how the prognostic FOX markers FOXD4, FOXH1, and FOXS1 correlated with the immune scores, stromal scores, tumor purity, and TIL quantities in the TCGA COAD.First, we found that the stromal scores, immune scores, and tumor purity are significantly deregulated (Wilcoxon sum rank test, P < 0.05) between the HEG and LEG of FOXH1 and FOXS1 (Figure 6A).However, the expression of FOXD4 is not significantly associated with these immune factors (Wilcoxon sum rank test, P < 0.05) (Figure 6A).Second, we investigated the correlation (Spearman correlation test, R > 0.20, P < 0.01) of these three independent  6B).Our analysis revealed that FOXS1 level showed a positive relationship with different immune signatures, such as CAFs, endothelial cells, HLA genes, immune cell infiltrate genes, immune checkpoint genes, macrophages, M2 macrophages, MDSCs, neutrophils, pDCs, T cell exhaustion, TAM, Tfh, TILs, Treg, type I IFN response, and Type II IFN response (Spearman correlation test, R > 0.20, P < 0.01) (Figure 6B).The FOXH1 level showed a negative correlation with the CD8 + regulatory T cells score, cytolytic activity, NK cells, and T cell activation (Spearman correlation test, R > 0.20, P < 0.01).However, the FOXD4 level was not related to the tumor immunity in the COAD (Spearman correlation test, R > 0.20, P < 0.01) (Figure 6B).Third, based on the level of FOXH1 and FOXS1 genes, we found that the ratios of pro-/antiinflammatory cytokines were lower in the HEG patients compared to the LEG patients (Wilcoxon sum rank test, P < 0.005) (Figure 6C).Fourth, we discovered a positive correlation between different immune-inhibitory hallmark markers and the level of FOXS1 expression in COAD (Pearson's correlation test, R > 0.20, P < 0.05) (Table 2).Some of the immune-inhibitory prominent marks, including IL10, PD-1, PD-L1, PD-L2, CTLA4, HAVCR2, TIGIT, CXCL13, FAP, VIM, and POSTN, were positively correlated with the expression level of FOXS1 (Table 2).Altogether, our immunological investigations demonstrated that the overexpressed level of FOXH1 and FOXS1 is crucially related to the immunosuppression of the colon TME, which may provide a significant explanation for the oncogenic role of FOX genes.

Prognostic FOX Genes Associated with Cancerassociated Pathways and Biological Processes in COAD
After we identified that the expression of three FOX genes (FOXD4, FOXH1, and FOXS1) are independent prognostic factors and regulate tumor immunity in COAD, we investigated the Spearman's correlation of these three genes with various cancer-associated pathways activity (R > 0.20 and P < 0.01).Interestingly, we discovered a favorable correlation between the level of FOXS1 expression and the activities of several pathways, such as colorectal cancer, ECM receptor interaction, ERBB, Focal adhesion, Hippo, JAK-STAT, MAPK, mTOR, Notch, Pathways in cancer, TGF beta, VEGF, WNT, and Hedgehog signaling pathway (Figure 7A).In contrast, the level of FOXS1 in COAD is inversely linked with cell cycle activity (Figure 7A).Moreover, we studied the association between the expression of these three distinct prognostic variables and biological processes associated with cancer, such as EMT, angiogenesis, and hypoxia.We revealed a favorable correlation between EMT, angiogenesis, hypoxia, and the level of FOXS1 in COAD (Figure 7B).This implies that FOXS1 expression crucially controls the biological processes and mechanisms involved with malignancy in COAD.

Genetic Alteration of Prognostic FOX Genes in Colorectal Cancer
To identify the genetic changes of three prognostic FOX genes (FOXD4, FOXH1, and FOXS1), we utilized the COAD (TCGA, Firehose Legacy) datasets with mutation and CNA information (n = 220) in cBioPortal (http:// www.cbioportal.org/).The 37 (17%) patient populations with the studied FOXD4, FOXH1, and FOXS1 genes showed genetic alterations.We revealed that FOXD4 was genetically changed in 2.7% of samples, and FOXH1 was genetically changed in 5% of samples (Figure 8A).In addition, FOXS1 was altered in 10% of patients (Figure 8A).The genetic changes mainly included amplification and missense mutation of three prognostic FOX genes (Figure 8A).

Expression Evaluation and Diagnostic Efficacy of Prognostic FOX Genes in Colorectal Cancer
We examined the mRNA expression of three prognostic FOX genes (FOXD4, FOXH1, and FOXS1) in a variety of malignancies, including colorectal cancer, using the Oncomine repository (https://www.oncomine.org/resource/login.html)(Figure 1 and Table 1).The value of FOXD4, FOXH1, and FOXS1 were significantly upregulated in colorectal cancer relative to the control (FC > 1.5 and P < 0.05) (Figure 8B).To validate the expression levels of these three genes, we used GSE166427 15 to differentiate the expression level of FOXD4, FOXH1, and FOXS1 between the COAD and normal adjacent colon cells.Interestingly, we found that FOXD4, FOXH1, and FOXS1 are consistently deregulated in GSE166427 15 (Figure 8C).
Since these three prognostic FOX gene members (FOXD4, FOXH1, and FOXS1) are key regulatory genes in COAD, we hypothesized that the three genes (FOXD4, FOXH1, and FOXS1) are useful for diagnosing COAD patients.We tested our stated hypothesis using the TCGA COAD dataset.For TCGA-COAD and healthy samples, the ROC curve of the expression levels The ratios of pro-/ anti-inflammatory cytokines were lower in the HEG patients compared to the LEG patients (Wilcoxon sum rank test, P < 0.005).We used the IFNG, IL-1A, IL-1B, and IL-2 marker genes as the pro-inflammatory cytokines and the IL-4, IL-10, IL-11, and TGFB1 marker genes as the anti-inflammatory cytokines 27 .* is P < 0.05, ** is P < 0.01, and *** is P < 0.001 of FOXD4 (AUC = 0.779), FOXH1 (AUC = 0.77), and FOXS1 (AUC = 0.861) displayed remarkable diagnostic significance (Figure 8D).Also, the ROC curve using the GSE166427 showed consistent diagnostic efficacy (Figure 8E).

Discussion
Colorectal cancer is the third leading diagnosed cancer with 9.4% of deaths worldwide. 143][4][5] Here, we explored the aberrant levels of FOX genes in COAD, the involvement of deregulated FOX genes in the functional enrichment, and the association of FOX genes in the PPI.Moreover, the FOX genes are correlated with poor clinical outcomes, tumor immunity, cancer-associated signaling pathways, and cancer hallmark biological processes.The COAD's 32 FOX genes were deregulated (Table 1 and Figure 1).Based on earlier research, colorectal tumors frequently display different unregulated FOX genes.For example, FOXQ1, the top up-regulated gene with Log2FC 5.86, is upregulated in colorectal cancer tissue samples and promotes cancer metastasis by regulating the PI3K/AKT signaling pathways. 36The elevated expression level of FOXD1 was found in human CRC tissues, and FOXD1 expression levels were correlated with tumor size and other clinical factors. 37FOXA2, another significantly upregulated gene, is overexpressed in colon cancer tissues, promotes EMT, inhibits apoptosis, and enhances the invasion ability of colon cancer cells. 38FOXP2, the top down-regulated FOX genes in COAD, promotes the invasion of hepatocellular carcinoma. 39Through controlling the EGFR/RAS/RAF/ MEK/ERK axis, the second-ranked down-regulated gene, FOXD3, serves as a tumor suppressor to reduce the malignancy in colon tissue. 40The expression of FOXF2 in Hela cells and tumor tissues was lower than nearby tissues, and FOXF2 has been linked to preventing Hela cells from proliferating, migrating, and invading through controlling the Wnt biochemical signaling process. 41onsistent with these studies, our findings revealed that the deregulated FOX genes are crucially associated with colon carcinogenesis.
Since the PPI network is considered an anticancer therapeutic discovery, 42 we built the PPI network of deregulated FOX genes.Interestingly, we found the five sub-networks (Figure 3) associated with FOX genes.Subnetwork 1 (Figure 3A) is crucially linked to the various cancer-associated signaling processes (Figure 3A).Numerous pathways consistently regulate colorectal cancer and other malignancies.Farhan et al. reported that the FOXO signaling pathways are substantial therapeutic targets in cancer treatment strategy. 43Uddin et al. established that the various pathways in cancer, small cell lung cancer, p53, apoptosis, and notch signaling processes are enriched in the colon tumor microenvironment. 44OXO subfamily members regulate the PI3K-AKT molecular pathway. 7The disruption of TGF-β signaling is crucially correlated with CRC tumorigenesis and progression. 45Altogether, the FOX gene-associated PPI network-mediated signaling cascades may be associated with colorectal carcinogenesis.
We identified that shorter COAD patient survival time is related to eight FOX genes (Figure 4).In numerous malignancies, the patient's survival duration consistently correlates with the prognostic FOX genes we uncovered.
For example, worse OS and DFS in colorectal cancer are linked to FOXC1 overexpression. 46Individuals with colorectal cancer have a poor prognosis when FOXD1 is up-regulated. 37Chen et al found that FOXD4 activity was increased in CRC and correlated with a short survival period. 47Reduced lung cancer survival rates are linked to up-regulated FOXH1 levels. 48FOXQ1, another prognostic gene, is poorly associated with prognosis in NSCLC. 49In patients with gastric cancer, overexpression of FOXS1 is related to a shorter survival period. 50Moreover, univariate, multivariable, and nomogram analyses showed that FOXD4, FOXH1, and FOXS1 are independent prognostic markers in COAD (Figure 5).These results indicate that those with colorectal cancer have a significantly reduced survival probability due to the dysfunctional activity of FOX genes.
The survival time of cancer patients is connected with immunogenicity, and immunological responses have a considerable effect on the clinical outcomes of patients. 51e evaluated the correlation of three independent prognostic FOX genes (FOXD4, FOXH1, and FOXS1) with tumor immunity in COAD (Figure 6).We found that the FOXH1 and FOXS1 genes substantially regulate tumor immunity in COAD (Figure 6).First, the immune and stromal content are deregulated between the HEG and LEG of these FOX genes (Figure 6A).Second, these two prognostic FOX genes are significantly correlated with immune infiltrations (Figure 6B).Third, the FOXH1 and FOXS1 highly expressed group in COAD had considerably reduced ratios of pro-/anti-inflammatory cytokines (Figure 6C).Fourth, the FOXH1 and FOXS1 genes are correlated with immunosuppressive markers (Table 2).Tumor-associated macrophages (TAMs), leading immunosuppressive immune cells, stimulate CRC growth by modifying the extracellular matrix remodeling, tumor metabolism, angiogenesis, and the tumor microenvironment. 52In metastasis of colorectal carcinoma, M2 macrophages are associated with carcinogenesis and tumor development by enhancing the invasiveness of tumor cells. 53CAFs, another immunosuppressive member, 54 improve the enrichment for TAMs and suppress NK activity in colorectal cancer. 55he immunosuppressive MDSCs are crucially related to the initiation and progression of CRC. 56The regulatory T cells (Tregs) contribute to the progression of human colorectal cancer. 57T cell exhaustion biomarkers, such as TIM-3/ HAVCR2, PD-1, CTLA4, and TIGIT, were overexpressed in CRC tumor tissues. 58Altogether, the prognostic FOX genes may be associated with the immunosuppressive colon tumor microenvironment by regulating the immune-inhibitory markers, stromal components, and immune cells.
We also evaluated the relationships between prognostic FOX genes and pathways linked to cancer and biological processes (Figure 7).Previous studies consistently found that the correlated cancer-associated pathways are enriched in colorectal cancer.For example, Uddin et al. found that cancer-related signaling processes, such as ECM-receptor interaction, focal adhesion, and Wnt signaling pathway, are significantly enriched in the colon tumor microenvironment. 44Hedgehog signaling is associated with CRC tissue formation, proliferation, and metastasis. 591][62] These suggest that the deregulated prognostic FOX genes are associated with COAD carcinogenesis by regulating the cancer-associated pathways and hallmark biological processes.
Ultimately, we evaluated the genetic alteration and diagnostic efficacy of three independent prognostic factors (FOXD4, FOXH1, and FOXS1) in COAD (Figure 8).In addition, we investigated the expression variation of these prognostic markers in cancers using the Oncomine database and GSE166427 (Figure 8), and found consistent results which ultimately indicate the substantial roles of FOXD4, FOXH1, and FOXS1 in colon cancer.Furthermore, molecular docking studies revealed that FOXH1 potentially interacted with various drugs (Figure 9 and Figure 10).Li et al reported that FOXD4 was consistently found to have diagnostic and prognostic value in colonic cancer and to be linked to the nuclear matrix, Rap1 signaling pathway, RNA transportation, and VEGF signaling framework. 63The growth of gastric cancer is connected with aberrantly expressed FOXS1, which has the potential to serve as a biomarker for both diagnosis and prognosis. 50Overall, these findings indicate that these FOX genes may be crucial targets for the therapeutic implications of COAD.
Our study also has some drawbacks.First, our findings were obtained from a bioinformatics study and need experimental validation.Second, we analyzed COAD's mRNA expression profiles, which are not equivalent to the protein expression levels.Hence, additional clinical and experimental verification will be required to convert these discoveries into practical uses for therapeutic targeting.
In conclusion, the dysfunctional FOX genes are substantially involved with poor clinical outcomes, immunogenicity, cancer-associated pathways, and oncogenic biological processes in COAD.This study may provide novel clues for targeting the FOX genes as potential therapeutics in COAD.

Figure 3 .
Figure 3. Involvement of FOX Genes in the PPI Network and the Relationship of FOX Genes with Other Family Members. A. Involvement of FOX genes in the PPI (sub-network 1) network.B. The top 30 KEGG pathways are considerably enriched and connected to sub-network 1. C-F.The other four sub-networks of FOX genes in the PPI.G. Relationship between FOX genes in the TCGA COAD data (Pearson correlation test).× represents the non-significant value (P > 0.05)

Figure 4 .
Figure 4. Exploration of Prognostic FOX Genes in the COAD. A. The HEG patients of FOXC1, FOXD1, FOXD4, FOXD4L2, FOXH1, FOXQ1, and FOXS1 showed a substantially poor OS rate in COAD.The LEG patients of FOXN1 showed a substantially poor OS rate in COAD.B. The HEG patients of FOXC1 and FOXD4L2 showed a substantially shorter DFS rate in COAD

Figure 5 .
Figure 5. Univariate and Multivariable Analysis of OS for Identifying the Prognostic Factors Associated with Nomogram Prediction. A. The univariate Cox model explored the risk factors in COAD, including FOXC1, FOXD1, FOXD4, FOXH1, and FOXS1.B. The multivariable Cox found that the stage, fraction genome altered, FOXD4, FOXH1, and FOXS1 are risk factors in COAD.C. Identifying the risk factors in the constructed nomogram.The nomogram predicted COAD's 1, 2, and 3 years OS.Fraction = fraction genome alterations

Figure 6 .
Figure 6.Prognostic FOX Genes Correlated with the Regulation of Immunity in the COAD. A. Immune scores, stromal scores, and tumor purity are deregulated with the expression of prognostic FOX genes.B. The correlation of various immune cells with the prognostic FOX genes in the COAD.C.The ratios of pro-/ anti-inflammatory cytokines were lower in the HEG patients compared to the LEG patients (Wilcoxon sum rank test, P < 0.005).We used the IFNG, IL-1A, IL-1B, and IL-2 marker genes as the pro-inflammatory cytokines and the IL-4, IL-10, IL-11, and TGFB1 marker genes as the anti-inflammatory cytokines27 .* is P < 0.05, ** is P < 0.01, and *** is P < 0.001

Figure 7 .
Figure 7. Association of Biological Signaling with the Expression of Prognostic FOX Genes in COAD. A. The cancer-associated pathways are correlated with prognostic FOX genes.B. The hallmark biological processes, including EMT, angiogenesis, and hypoxia, correlate with the expression level of FOXS1.* is P < 0.05, ** is P < 0.01, and *** is P < 0.001

Figure 8 .
Figure 8. Genetic Alterations, Expression Validation, and Diagnostic Efficacy of Prognostic FOX Genes in Colorectal Cancer. A. Genetical alterations of FOXD4, FOXH1, and FOXS1 in TCGA COAD cohort.The FOXD4 (2.7%), FOXH1 (5%), and FOXS1 (10 %) genes are mutated in COAD.B. mRNA expression of the FOXD4, FOXH1, and FOXS1 genes in different malignancies (results retrieved from the Oncomine repository).The fold change was set to 1.5, and P-value was set at 0.05.The value in the table indicates how many datasets meet the criteria.The extent of overexpression or decreased expression is linked with the strength of the color (red or blue).A light-green box indicated the mRNA deregulation of three prognostic genes.C. Consistent expression of FOXD4, FOXH1, and FOXS1 in GSE166427.D. Diagnostic efficacy of three prognostic FOX genes in TCGA-COAD.We sketched the receiver operating characteristic (ROC) curve of three prognostic genes in COAD and healthy samples.The value of FOXD4, FOXH1, and FOXS1 revealed diagnostic efficacy in the COAD.E. Consistent diagnostic efficacy of three prognostic FOX genes in GSE166427 15

Table 1 .
Differential Expression of FOX Genes in TCGA COAD

Table 2 .
The Independent Prognostic FOX Genes Correlated with Immunosuppressive Markers in the COAD R is Pearson's correlation, and P is significance level.