Expression of the HOXA gene family and its relationship to prognosis and immune infiltrates in cervical cancer

Abstract Background The homeobox A cluster (HOXA) gene family is participated in multiple biological functions in human cancers. To date, little is known about the expression profile and clinical significance of HOXA genes in cervical cancer. Methods We downloaded RNASeq data of cervical cancer from The Cancer Genome Atlas (TCGA) database. The difference in HOXA family expression was analyzed using independent samples t test. Cox proportional hazard regression analysis was used to assess the effect of HOXA family expression on survival, and a nomogram predicting survival was generated. We assessed the infiltration difference in immune cells and expression difference of immunity biomarkers between two groups with different expression level of HOXA genes through Immune Cell Abundance Identifier (ImmuCellAI) and independent samples t test, respectively. Results Our results showed that the HOXA1 gene was upregulated, while the HOXA10 and HOXA11 were downregulated in cervical cancer. Downregulation of HOXA1 was related to a poor outcome for cervical cancer patient. We also identified a significantly increased abundance of T helper 2 cells (Th2) and higher expression of PD‐L1 in cervical cancer patients with lower expression of HOXA10 and HOXA11. The gene set enrichment analysis (GSEA) results indicated that HOXA1 and HOXA11 were involved in immune responses pathways and participated in the activation of a variety of classic signaling pathways related to the progression of human cancer. Conclusion This study comprehensively analyzed different HOXA genes applying public database to determine their expression patterns, potential diagnostic, prognostic, and treatment values in cervical cancer.

of human papillomavirus (HPV) for cervical cancer, H. pylori for stomach cancer, and HBV for liver cancer. The proven biological etiology of cervical cancer is HPV infection, in which HPV 16 and 18 are attributed to more than 70% of cases worldwide. 3 The brush over of sexual education and economic backwardness of China are responsible for HPV infection. The most effective strategy for the prevention of infection-related cancers is to create more effective vaccines against these carcinogenic viruses and to formulate better annihilation methods to combat these bacteria. In addition, appropriate screening methods are critical to improve the overall survival (OS) of cases. The adherence rates in the cervical cancer screening program were relatively high in the USA, with an adherence rate of 83%. 4 At present, surgery, concurrent chemoradiotherapy and radiotherapy have better effects on early cervical cancer. Most of early stage of cervical cancer (stage I to IIa) can be cured by surgery or radiotherapy, while treatment should be tailored for advanced and metastatic cervical cancer, and concurrent chemoradiation is considered as the priority treatment. However, the efficacy is still limited. 5 Therefore, the identification of biomarkers for early cervical cancer could be propitious to favor the outcome of cervical cancer.
Biomarkers play crucial roles in predicting the treatment response, prognosis, and disease progression in cancer, developing new therapies, and elucidating tumorigenesis mechanisms. 6 The homeobox genes (HOX) are a series of genes coded crucial transcriptional regulators that play diverse roles from embryogenesis to tumorigenesis. 7 In humans, the HOX family is arranged into four gene clusters termed HOXA, HOXB, HOXC, and HOXD. 7 HOXA cluster genes include HOXA1~7, HOXA9~11, and HOXA13. 7 Increasing evidence has shown that altered HOXA genes are involved in human cancer initiation and progression. HOXA13 was expressed more in the normal colon tissues when compared with colon cancer, and HOXA13 expression was associated with TNM stage of patients. 8 HOXA13 is also involved in HOTTIP-induced malignant phenotypes of gastric cancer cells. 9 It was also reported that HOXA2 and HOXA4 were downregulated and other HOXA members were upregulated in laryngeal squamous cell cancer, including HOXA9 and HOXA13. 10 In addition, upregulation of HOXA10, HOXA11, and HOXA13 was correlated with poor laryngeal squamous cell cancer OS. 10 However, to date, little is known about the expression profile and clinical significance of HOXA genes in cervical cancer.
In the present study, we analyzed different HOXA genes applying public database to determine their expression patterns, diagnostic, prognostic potential, and treatment values in cervical cancer.

| The cancer genome atlas (TCGA) data acquisition
The Cancer Genome Atlas (TCGA) projects (https://portal.gdc.cancer. gov/), consisting of excellent databases of very large RNA sequence data of cancerous and normal samples, provide great opportunities for high-throughput modeling and bioinformatics analysis to determine diagnostic and prognostic biomarkers of cancer. The TCGA cervical cancer database includes a total of 306 cervical cancer and 3 normal control samples. We obtained high-throughput sequencing (HTSeq) data, Illumina Human Methylation 450K data, and clinical parameters (including TNM stage, smoking history, grade, and age).

| Expression and methylation data of HOXA genes in cervical cancer
The expression and methylation data of the HOXA genes were abstracted from HTSeq data and Illumina Human Methylation 450K data. The comparison of expression level of all HOXA genes was analyzed using an independent samples t test between cervical cancer and normal tissues.

| Receiver operating characteristic (ROC) curve analysis
We performed a receiver operating curve (ROC) analysis to evaluate the diagnostic value of determination of HOXA expression in distinguishing cervical cancer patient and to calculate the area under the curve (AUC) value and cutoff value based on the maximum of Youden index. 11 The Youden index is defined as sensitivity +specificity −1. 11 Then, cervical cancer patients were grouped into two groups according to the cutoff of HOXA expression, respectively. Then, time-dependent receiver operating characteristic (ROC) curve analyses and AUC were also performed to evaluate the prediction accuracy of our model.

| Prognostic values of HOXA members in cervical cancer
The prognostic potential of the HOXA genes was assessed by Cox proportional hazards regression analyses. First, the associations between all clinical parameters (including HOXA genes expression, TNM stage, age, smoking history, and histological grade) and overall survival (OS) among cervical cancer patients were assessed using univariate Cox proportional hazards regression analyses.
Subsequently, all the survival associated clinical features, including HOXA expression, were evaluated using multivariate Cox proportional hazards regression analysis.

| Construction of nomograms
A nomogram generates a numerical probability of a clinical event based on a series of clinical features and a statistical predictive model. 12 The development of a nomogram includes defining the patient outcomes, identifying important covariates, specifying the statistical model, and validating its performance. Briefly, the outcome is typically an event, which is the time to death in the current study.
Nomograms were used to predict the probability of the 3-and 5- year OS. The final model selection was performed using a backward step-down process in the Cox proportional hazards regression analysis, that is, the prognostic covariates included the significant prognostic parameters assessed in our study.
Subsequently, a time-dependent receiver operating characteristic (ROC) curve analysis involving the 3-and 5-year OS was conducted to calculate the prognostic accuracy of the model for time-dependent survival. The performance of the nomogram was measured using a concordance index (C-index), which quantifies the level of concordance between the predicted and actual survival probabilities. The calibration was evaluated by plotting the relationship between actual and predicted probabilities using a bootstrapping method. 13

| Negative association between expression and methylation of the HOXA genes in cervical cancer
We also obtained the methylation levels of cg sites in the gene promoter regions of differentially expressed HOXA genes in cervical cancer. Then, we performed the Pearson's analysis to assess the association between HOXA expression and methylation in cervical cancer and applied the corrplot package to plot the association.

| Immune cell abundance analysis
Immune Cell Abundance Identifier (ImmuCellAI) is an online tool to estimate the abundance of 24 immune cells from a gene expression dataset including RNA-Seq and microarray data, in which the 24 immune cells comprise 18 T-cell subtypes and 6 other immune cells: B cells, NK cells, monocytes, macrophages, neutrophils, and DCs. 14 We used this tool to assess the infiltration difference in immune cells between the low-and high-HOXA expression groups of cervical cancer. Besides, we also compared the expression of relevant immunity biomarkers (including programmed cell death protein 1 (PD-1), programmed cell death ligand 1 (PD-L1), and cytotoxic T lymphocyte-associated antigen-4 (CTLA4)) between low and high expression of HOXA groups.

| Gene set enrichment analysis (GSEA)
The and HOXA11 groups were performed. Pathways with false discovery rate (FDR) <0.05 were considered with statistical significance.

| Statistical analysis
The HTSeq data and methylation data of HOXA genes from the TCGA database were abstracted using Perl 5.28. The association between methylation and expression of HOXA members were assessed and plotted by the corrplot package, the survival package was used for the analysis of prognostic values, and the rms package was applied for the construction of the nomogram. The comparison of expression of HOXA genes in cervical cancer tissues and normal controls was analyzed using independent samples t test by SPSS 25.0.

| Expression profile of HOXA family in cervical cancer
To obtain a full picture of the expression status of HOXA members in cervical cancer, we downloaded the HTSeq data of HOXA expression, which originated from the TCGA database (https://portal.gdc. the HOXA members (including HOXA2-9 and HOXA13) showed very small differences between cervical cancer and normal controls, with no statistical significance ( Figure S1).

| Diagnostic value of HOXA members in cervical cancer
To evaluate the discriminative ability of HOXA members, a receiver

| Prognostic value of HOXA members in cervical cancer
Furthermore, the prognostic values of HOXA members were estimated. Firstly, the survival package of R software was utilized to execute the Kaplan-Meier plot. The results were exhibited in Figure 3 and revealed that low expression of seven HOXA members were correlated with favor overall survival of cervical cancer (including HOXA1, HOXA2, HOAX3, HOXA4, HOXA5, HOXA6, and HOXA9).
Considering the principle of the survival package (seeking the optimal cutoff for overall survival, grouping patient in two group, and then drawing the Kaplan-Meier plot), we then determined the pre-

| Prognostic nomogram for OS in cervical cancer
Predictive models with nomograms were constructed using the rms package in R software, integrating age, histological grade, TNM stage, and HOXA expression, as shown in Figure 4. Each prognostic parameter was assigned a score according to its prognostic value; the sum total of the scores was used to predict 3-and 5-year OS. The total score for all the variables was converted into an estimate of the probability of death. The performance of the nomogram was measured by the concordance index (C-index); the larger the C-index, the more accurate the prognosis. The C-index for overall survival pre-

| Negative correlation of HOXA expression and promoter methylation
To create multiple views of the roles of the HOXA family in cervical cancer, we simultaneously obtained the methylation data of HOXA members. We identified three differentially expressed HOXA

| Immune cell abundance analysis
Firstly, patients with cervical cancer were classified into two different group based on the cutoff value of HOXA1, HOXA10, and HOXA11 expression from ROC. We analyzed the infiltration difference in immune cells between the low-and high-HOXA1, HOXA10, and HOXA11 expression groups using the ImmuCellAI online tool. We found that the abundance of T helper 2 cells (Th2) was significantly increased in both the low-HOXA10 and HOXA11 expression groups, while the abundance of dendritic cells (DCs) was increased in the low-HOXA11 expression group ( Figure 7B,C). There was no difference in the infiltration of immune cells between the low-and high-HOXA1 expression groups ( Figure 7A). Subsequently, we also analyzed the difference expression of several immunity biomarkers (CTLA4, PD-1, and PD-L1) between the low-and high-HOXA1, HOXA10, and HOXA11 expression groups. Our results showed that only patients with low expression of HOXA10 and HOXA11 were characterized with relatively high expression of PD-L1 (Figure 8).

| Analyses of the signaling pathways participated in cervical cancer using GSEA
The GSEA results revealed that multiple cancer-related pathways (including NOTCH signaling, p53 signaling, apoptosis pathway, and pathways in cancer) were enriched in the cervical cancer patients with high expression of HOXA1 and low expression of HOXA11 ( Figure 9, detail information shown in Tables S1 and S3). Among  (Table S2).

| DISCUSS ION
Homeobox (HOX) genes were firstly observed in the fruit fly, a large which are located on various chromosomes. 16 Each cluster is divided into 13 paralog groups (HOXA1-13 genes) based on the position of chromosomes. 17 These genes are characterized by a consensus DNA sequence. 18 Evidence of multiple studies has identified biological functions for HOXA genes during morphogenesis, patterning, and differentiation. [19][20][21] It was reported that HOXA is the most highly expressed in the developing lungs of mice 22,23 ; however, mice with mutations in HOXA5 at birth lack pulmonary surfactant, indicating a role of HOXA5 gene in morphogenesis. 24 Abnormal expression of HOXA genes has been proven in multiple cancers. 25,26 The oncogenic potential of HOXA genes has been clearly implicated in leukemia, and their roles in other neoplasia are also being evaluated. 27 The increased expression of HOXA9 could be found in most aggressive leukemia, 28 and HOXA13 physically linked to tumor growth and angiogenesis. 29 Numerous reports have cataloged differences in HOXA gene expression between normal and neoplastic tissues. HOXA1 and HOXA7 were expressed in small cell lung cancer, while they were not expressed in normal lungs.
Therefore, the authors proposed that the HOXA gene profiles of cells obtained from sputum could act as molecular markers to aid the detection of lung carcinoma. 30  cancer. In addition, our results showed that the expression level of HOXA members (HOXA1, HOXA10, and HOXA11) are affected by the methylation level, which has been reported by previous study. 10,32 A previous study identified HOXA4 as the only prognostic gene in acute myeloid leukemia (AML), the upregulation of which in patients with normal karyotypes was actually associated with longerterm survival. 28 HOXA13 may serve as a prognostic parameter in kidney renal clear cell carcinoma patients. 33 In our research, we also assessed the prognostic ability of the HOXA family in cervical cancer  (HOXA1, HOXA10, and HOXA11), HOXA1 was the only prognostic gene in cervical cancer. However, among the remaining members of the HOXA family, HOXA3 also exhibited prognostic capability without differential expression in cervical cancer. We also developed a clinical tool (a nomogram) that predicts  36,37 In this study, we applied the online tool ImmuCellAI, a highly accurate method of estimating the abundance of immune cells. 14

CO N FLI C T O F I NTE R E S T
The authors declare that there are no conflicts of interest.

AUTH O R CO NTR I B UTI O N S
FG and JZ were responsible for designing the research study. YZ, WT, and YF were in charge of analyzing the data from public database. FG was responsible for writing of article. All authors reviewed the study.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in TCGA database at https://portal.gdc.cancer.gov/.

F I G U R E 9
The Gene Set Enrichment Analysis (GSEA) of the relationship between the expression level of HOXA1 and HOXA11 in cervical cancer