A total of 53 NSCLC patients, including 31 lung adenocarcinoma (LUAD) patients and 22 lung squamous cell carcinoma (LUSC) patients, were enrolled from the Department of Pathology, First Affiliated Hospital of Guangxi Medical University (Nanning, Guangxi, China). All 53 samples were randomly collected from patients undergoing surgical resection without treatment. All methods were applied according to the relevant guidelines. The Ethics Committee of the First Affiliated Hospital of Guangxi Medical University approved the experimental protocols, and all patients provided written informed consent forms for the use of their tissues in this study. Total RNA was extracted via TRIzol reagent (Thermo Fisher Scientific, Inc., Waltham, MA, USA), and a PCR amplification kit (Omega, Solarbio Biotechnologies, Inc., Shanghai, China) was used. The RNA was reverse-transcribed into cDNA using Roche cDNA Synthesis kit (Roche Diagnostics, Shanghai, China), based on the manufacturer's protocols. qPCR was performed using ABI 7500 prepstation (Applied Biosystems; Thermo Fisher Scientific, Inc.), and the SYBR^®-Green PCR Master mix (GeneCore Biotechnologies, Inc., Shanghai, China). PCR was performed at 95°C for 15 sec, 60°C for 1 min, 95°C for 15 sec and 60°C for 1 min for 40 cycles. The specific primers were as follows: HOXA1 forward, 5′-CGGCTTCCTGTGCTAAGTCT-3′ and reverse, 5′-TAGCCCAGCCAAATACACGG-3′; and GAPDH (internal control) forward, 5-′TGCACCACCAACTGCTTA-3′ and reverse, 5′-GGATGCAGGGATGATGTTC-3′. The results were normalized to the GAPDH expression and calculated based on the 2^−ΔΔCq method (17,18).

TCGA (http://cancergenome.nih.gov/) has collected comprehensive molecular profiles, including gene expression, microRNA expression, protein expression and DNA methylation, for >30 types of human tumors (19–21). TCGA also has information about complex clinical parameters. In the present study, the RNA-Seq data for patients with NSCLC, which were from the Illumina HiSeq RNA-Seq platform (Illumina, Inc., San Diego, CA, USA), contained 535 LUAD cases and 502 LUSC cases up to July 1, 2017 (21). The expression data of HOXA1 are reported in reads per million, and the HOXA1 expression level was normalized by the R language package DESeq for further analysis. Student's t-test (SPSS Inc., Chicago, IL, USA) was used to compare differential expression of HOXA1 between NSCLC and normal lung tissues. Additionally, the potential associations between HOXA1 and the clinicopathological parameters in NSCLC were identified via the original TCGA database. The receiver operating characteristic (ROC) curve was derived to evaluate the diagnostic value of HOXA1. Oncomine (https://www.oncomine.org/) and GEPIA (http://gepia.cancer-pku.cn/) were applied to verify the HOXA1 expression in NSCLC (22,23).

To further investigate the genes co-expressed with HOXA1, MEM (http://biit.cs.ut.ee/mem/index.cgi), GEPIA and cBioPortal (http://www.cbioportal.org/) were used. The Venn diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/) were used to identify and compare the overlaps. Next, bioinformatic analyses, including GO, KEGG and network analyses, were utilized to investigate the potential functions, pathways and networks of these overlapping genes as previously described (24). In this process, the Database for Annotation, Visualization and Integrated Discovery (http://david.abcc.ncifcrf.gov/) was used for GO and KEGG analyses. Biological process, cellular component and molecular function were derived separately via GO analysis. A functional network was constructed through Cytoscape (version 2.8; http://cytoscape.org).

The interaction pairs of the co-expressed genes were researched through the Search Tool for the Retrieval of Interacting Genes (STRING; version 9.0; http://string-db.org) (25). The STRING database aims to supply a global perspective for as many organisms as feasible. Known and predicted associations are integrated and scored. A combined score over 0.4 was chosen to construct the PPI network.

A total of 12 target prediction algorithms were used for predicting the potential miRNAs targeting HOXA1: miRWalk (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/), DIANA microT v4 (http://diana.imis.athena-innovation.gr/), miRanda (http://www.microrna.org), mirBridge (http://mirsystem.cgm.ntu.edu.tw/), miRDB (http://www.mirdb.org/), miRMap (http://mirmap.ezlab.org/), miRNAMap (http://mirnamap.mbc.nctu.edu.tw/), Pictar2 (https://www.mdc-berlin.de/), PITA (https://genie.weizmann.ac.il/), RNA22 (https://cm.jefferson.edu/) RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/) and TargetScan (http://www.targetscan.org/). Candidate miRNAs were identified based on Venn diagrams.

All the original data from TCGA were log2-transformed. The mean ± standard deviation was calculated by SPSS 22.0 (IBM Corp., Armonk, NY, USA) to measure the HOXA1 expression level. Student's t-test was used to compare the differential expression of HOXA1 between NSCLC and normal lung tissues, as well as for the associations between HOXA1 expression and the clinicopathological parameters. One-way analysis of variance was applied to compare different subgroups. The Mann-Whitney U test or Kruskal-Wallis H test was utilized for non-normally distributed variables. The associations between HOXA1 expression and miRNA expression were assessed by Spearman's correlation. Mining for co-expressed genes across hundreds of datasets was performed through novel rank aggregation and visualization methods. Two-sided P-values of <0.05 were identified to indicate statistical significance.

In the present study, HOXA1 mRNA was overexpressed in NSCLC compared with that in normal lung tissues (P=0.011; Fig. 1A). The associations between the expression of HOXA1 and different clinicopathological parameters were further investigated. HOXA1 expression was positively associated with advanced (stage III–IV) TNM, Tumor-Node-Metastasis (TNM) classification of malignant tumors (26) and presence of lymph node metastasis (LNM) (both P<0.001; Fig. 1B and C). No significant association was found between HOXA1 mRNA expression and any other clinicopathological parameter, including sex, tumor size and vascular invasion (Table I). In addition, the diagnostic value of the HOXA1 level in NSCLC was assessed by ROC curve, and the area under the curve (AUC) of HOXA1 was 0.656 [95% confidence interval (CI), 0.552–0.761; P=0.005; Fig. 1D). The expression of HOXA1 was also compared between LUAD and LUSC. The results were similar to those of NSCLC: HOXA1 was upregulated in LUSC (P=0.012; Fig. 2A), and HOXA1 expression was positively associated with the presence of LNM and an advanced TNM stage (III-IV, both P<0.05) in LUAD and LUSC (Fig. 2B–2E; Tables II and III). A moderate diagnostic value of the HOXA1 level was also found in LUAD (0.625; 95% CI, 0.502–0.748; P=0.057; Fig. 2F), although this was not significant, and in LUSC (0.700; 95% CI, 0.566–0.834; P=0.007; Fig. 2G). Comparison of HOXA1 expression between LUAD and LUSC tissues showed higher expression of HOXA1 in LUSC than LUAD.

To further research the differential expression of HOXA1 between NSCLC and non-cancerous lung tissues, original patient data was obtained from TCGA. Two NSCLC cohorts, which comprised i) 535 LUAD cases and 59 normal lung cases, and ii) 502 LUSC cases and 49 normal lung cases, were extracted. As a result, increased expression of HOXA1 was observed in LUAD and LUSC compared with that in normal lung tissues (both P<0.05; Fig. 3A and B). Regarding the clinicopathological parameters, no statistical significance was reached based on the TCGA database (Tables IV and V). The AUC of HOXA1 was 0.548 (95% CI, 0.498–0.599; P=0.002) for LUAD and 0.957 (95% CI, 0.940–0.974; P<0.001) for LUSC based on TCGA, which indicated a high diagnostic value of the HOXA1 level in LUSC (Fig. 3C and D). With regard to overall survival, no statistical significance was determined; a trend was observed in which low HOXA1 expression was associated with an increased survival time (97.11±11.49 months) compared with high HOXA1 expression (75.15±9.68 months) (P=0.098; Fig. 3E) in LUAD, and the opposite trend was noted in LUSC (P=0.795; Fig. 3F), indicating that high HOXA1 expression may be associated with increased survival time of NSCLC patients.

A total of 11 datasets [Hou Lung, Wachi Lung, Beer Lung, Stearman Lung, Garber Lung, Landi Lung, Bhattacharjee Lung, Su Lung, Talbot Lung, Selamat Lung and Okayama Lung (22)] in Oncomine were used to validate the HOXA1 expression. Bhattacharjee Lung showed an opposite trend to all other datasets, as HOXA1 expression was downregulated compared with that in the normal lung. The results from the other 10 datasets were consistent with the present RT-qPCR and TCGA findings (Fig. 4A and B). GEPIA was used to further confirm the high expression of HOXA1 in LUAD and LUSC compared with that in the non-cancerous lung tissues (Fig. 4C and D).

Based on GEPIA, TCGA and MEM, 1,264 overlapping co-expressed genes were selected (Fig. 5) for GO and KEGG pathway analyses. The strongly enriched GO functional terms were 'dGTP metabolic process', 'network-forming collagen trimer' and 'centromeric DNA binding' (Fig. 6; Table VI). The KEGG pathway most strongly associated with the HOXA1 co-expressed genes was 'the p53 signaling pathway' (Table VII). Altogether, the GO and KEGG pathway analyses indicated that HOXA1 may be associated with the biological mechanism of NSCLC.

A PPI network was constructed via STRING online, and a total of 3,250 PPI pairs with a combined score of >0.4 were noted. The map of the PPI network that involved 908 PPI pairs was chosen for further analysis, and its connectivity degree was >30 (Fig. 7). Trifunctional purine biosynthetic protein adenosine-3 (GART; degree=71) had the highest degree and most interactions, according to the PPI network.

A total of 17 genes (ZMAT3, CYCS, CHEK1, CDK6, SFN, SESN3, CCNB1, CCNE1, TP53I3, CDKN2A, CCNB2, CCND2, SERPINB5, DDB2, PERP, IGFBP3 and GADD45A) associated with the p53 signaling pathway were flagged by KEGG pathway analysis, and 4 genes (CDKN2A, RAD51, CHEK1 and GART) had a degree of connectivity of >50 in the PPI network. The genes shared in common by these two lists were cyclin-dependent kinase inhibitor 2A (CDKN2A) and checkpoint kinase 1 (CHEK1). The expression of the two genes in the original TCGA data was investigated and it was found that each was highly expressed in LUAD and LUSC compared with that in normal lung tissues (both P<0.001; Fig. 8A–D). Based on these results, we hypothesized that HOXA1 serves a vital role in NSCLC by co-expressing with CDKN2A and CHEK1.

In the present study, 12 target prediction algorithms were used to predict the potential miRNAs that targeted HOXA1. The miRNAs predicted by >10 algorithms were selected as the final candidate miRNAs. A total of 3 miRNAs (miR-181b-5p, miR-28-5p, miR-181d-5p) targeting HOXA1 were predicted by the 10 algorithms. Based on TCGA, miR-181b, miR-28 and miR-181d levels were found to be significantly upregulated in LUAD compared with those in non-cancerous lung tissues (all P<0.05; Fig. 9A–C). miR-181b and miR-181d levels were found to be significantly downregulated in LUSC, whereas miR-28 level was found to exhibit no significant difference in LUSC and normal lung tissues (Fig. 9D–F). Furthermore, the correlation between HOXA1 and these three miRNAs in NSCLC was compared using Spearman-test based on TCGA, and it was found that HOXA1 mRNA level was inversely correlated with miR-181b and miR-181d in both LUAD and LUSC (both P<0.05; Table VIII). However, miR-28 level was inversely correlated with HOXA1 mRNA level in LUAD (r=−0.010, P=0.827) and positively correlated in LUSC (r=0.057, P=0.216) (Table VIII). Since miR-181b and miR-181d expression was downregulated in LUSC tissues and inversely correlated with HOXA1 expression, we speculate that HOXA1 may be the direct target of miR-181b-5p or miR-181d-5p in LUSC, and may serve a significant role in NSCLC by regulating various pathways, particularly the p53 signaling pathway. However, the detailed mechanism should be verified by functional experiments.