Diagnostic significance and carcinogenic mechanism of pan‐cancer gene POU5F1 in liver hepatocellular carcinoma

Abstract Background The prognostic and clinicopathological significance of POU Class 5 Homeobox 1 (POU5F1) among various cancers are disputable heretofore. The diagnostic value and functional mechanism of POU5F1 in liver hepatocellular carcinoma (LIHC) have not been studied thoroughly. Methods An integrative strategy of meta‐analysis, bioinformatics, and wet‐lab approach was used to explore the diagnostic and prognostic significance of POU5F1 in various types of tumors, especially in LIHC. Meta‐analysis was utilized to investigate the impact of POU5F1 on prognosis and clinicopathological parameters in various cancers. The expression level and diagnostic value of POU5F1 were assessed by qPCR in plasma collected from LIHC patients and controls. The correlation between POU5F1 and tumor infiltrating immune cells (TIICs) in LIHC was evaluated by CIBERSORT. Gene set enrichment analysis (GSEA) was performed based on TCGA. Hub genes and related pathways were identified on the basis of co‐expression genes of POU5F1. Results Elevated POU5F1 was associated with poor OS, DFS, RFS, and DSS in various cancers. POU5F1 was confirmed as an independent risk factor for LIHC and correlated with tumor occurrence, stage, and invasion depth. The combination of POU5F1 and AFP in plasma was with high diagnostic validity (AUC = 0.902, p < .001). Specifically, the level of POU5F1 was correlated with infiltrating levels of B cells, T cells, dendritic cells, and monocytes in LIHC. GSEA indicated that POU5F1 participated in multiple cancer‐related pathways and cell proliferation pathways. Moreover, CBX3, CCHCR1, and NFYC were filtered as the central hub genes of POU5F1. Conclusion Our study identified POU5F1 as a pan‐cancer gene that could not only be a prognostic and diagnostic biomarker in various cancers, especially in LIHC, but functionally carcinogenic in LIHC.


| INTRODUCTION
Cancer has become a key influencing factor of morbidity and mortality in both developed and developing countries. 1 There will be an escalating trend of death rates caused by cancers in the future due to deficient cognition in the pathological processes and regulatory mechanisms of cancers. 2 Although the prognoses of cancers have been ameliorated through various therapeutic methods, the prognostic outcomes are invariably unsatisfactory in multiple kinds of cancers. Among all kinds of cancers, liver hepatocellular carcinoma (LIHC) plays one of the major roles. According to the cancer statistics data from 2020, LIHC ranks sixth in mortality among all cancers. 3 As the main diagnostic and prognostic biomarker for LIHC, the sensitivity and specificity of α-fetoprotein (AFP) were ungratified in the early diagnosis of LIHC. Consequently, urgent requirements are raised to find novel biomarkers as potential diagnostic indicators and therapeutic targets of LIHC.
Many studies have certified that cancer stem cells (CSCs) are associated with aggression, metastasis, and recrudescence in various cancers. In addition, several CSC markers have been proven to contribute to the poor prognosis of cancers, [4][5][6][7] indicating the significance of CSC markers in the prognosis of malignancies. However, due to the complexity of the regulatory network in tumor pathologic processes, the prognostic significance of CSC markers is not fully understood. With more in-depth studies on CSC markers, some of these markers may become important targets in cancer diagnosis, therapy, and prognosis. POU Class 5 Homeobox 1 (POU5F1) is a transcription factor of the POU family that binds an octameric sequence motif to activate the expression of downstream genes. 8 POU5F1 has been identified as one of the most important CSC markers and participates in stemness maintenance in various tumors. 9,10 Published literatures have certified that increased POU5F1 was correlated with clinicopathological features and prognosis not only in LIHC, but also in bladder carcinoma, non-small cell lung carcinoma, and oral squamous cell cancer. [11][12][13][14] POU5F1 may serve as an essential predictive factor for multiple cancers in the near future.
Although many studies have been performed, the prognostic significance of POU5F1 in cancers remains controversial, and the functions of POU5F1 in the regulatory network of tumors are not fully recognized. Some studies have come to different or even totally opposite conclusions regarding the prognostic value of POU5F1 and the role of POU5F1 in tumor development. For instance, He et al. showed that elevated POU5F1 in esophageal squamous cell carcinoma symbolized poor survival outcomes. 15 However, Ge et al. found that high expression of POU5F1 was connected with longer survival in esophageal squamous cell carcinoma. 16 The prognostic value of POU5F1 in LIHC was not statistically significant according to Qian et al 17 ; but was prominent in studies performed by Huang et al. 18 These disputes have not been settled in a reasonable way and the value of POU5F1 in tumor prognosis is still ambiguous. Current studies on the role of POU5F1 in LIHC mainly used tissue samples, hindering the clinical application of POU5F1 as a diagnostic biomarker due to its invasiveness. Studies focusing on the POU5F1 status in plasma could facilitate its promotion.
In this study, we adopted an integrative strategy of meta-analysis, bioinformatics, and wet-lab approach to explore the diagnostic and prognostic significance of POU5F1 in various types of tumors, especially in LIHC. First, we performed a meta-analysis and trial sequential analysis (TSA) with a large sample size to evaluate the significance of POU5F1 for survival prognosis in various cancers. LIHC was selected as the major target when combining the meta-analysis results with the survival analysis results from TCGA datasets. Then, we validated the POU5F1 expression level in plasma and evaluated the diagnostic value of POU5F1 in LIHC. Furthermore, a protein-protein interaction (PPI) network was constructed based on the co-expression genes of POU5F1 and central hub genes were recognized. Finally, a cell signal transduction diagram was drawn to clarify the potential functional pathways of POU5F1 in LIHC.

| Literature search strategy
We comprehensively retrieved PubMed, Embase, Web of Science, and Cochrane Library to search studies published from 1 January 2000 to 1 June 2019 with language limitation in English and screened studies reporting prognosis and clinicopathological features in cancer patients with aberrant expression of POU5F1. To increase search sensitivity, we used a strategy involving both Medical Subject Heading terms and free-text words. The search strategy was segmented into three parts: "POU5F1 transcription factor or octamer transcription factor 4 or octamer transcription factor 3", "neoplasms or malignant neoplasms or carcinoma", and "prognosis or prognostic factors or survival". We also manually browsed the references of retrieved articles to recognize more eligible studies that might have been missed by the search strategy.

| Literature inclusion and exclusion criteria
Published articles that met the seven criteria were enrolled: (a) evaluated the association between POU5F1 expression and clinical prognosis or clinicopathological parameters of cancers; (b) provided hazard ratios (HRs) and 95% CI or survival curves of POU5F1 relevant outcomes; (c) cohort studies (follow-up duration longer than 24 months); (d) whole paper was written in English; (e) available fulltext articles; (f) research on humans; and (g) sample size of cancer patients was no less than 20. The exclusion criteria included the following: (a) absence of essential data, such as detection methods of POU5F1 expression, survival analysis data, and accurate prognosis indicators; and (b) reviews, case reports, letters, conference abstracts, animal trials, or duplicate publications.

| Literature data extraction and quality evaluation
Each process of our research was strictly in conformity to preferred reporting items for systematic reviews and meta-analyses (PRISMA) guidelines. 19 Important features of the eligible cohorts were recorded, including first author; published year; nation; sample size; tumor category; age and gender of the patients; detection method and cut-off value for POU5F1; follow-up period; study design; clinicopathological parameters; outcome of interest, including overall survival (OS), diseasefree survival (DFS), disease-specific survival (DSS), and recurrence-free survival (RFS). The Newcastle-Ottawa Scale (NOS) was utilized to appraise the quality of the included cohorts. According to the NOS criteria, a cohort was considered of high quality when the total score was no less than 7. 20

| Trial sequential analysis
With new studies constantly enrolled in the cumulative meta-analysis, type I and type II errors might increase due to repetitive tests of significance, fragmentary data, and ambiguous publication bias. 21 Trial sequential analysis (TSA) can overcome these obstacles and estimate a priori information size (APIS), which is considered as the minimal sample size required to draw a reliable conclusion. When the cumulative Z-curve fails to cross the conventional boundary (Z = 1.96), it indicates that the result is farfetched. If the Z-curve crosses the conventional boundary but does not reach the TSA boundary, meaning the trials show false positive results. If the Z-curve crosses both the conventional boundary and TSA boundary but not the APIS, it suggests that more researches are needed to support the conclusion. If the Z-curve crosses all three boundaries, a reliable conclusion has been certified. We implemented TSA by maintaining two-sided α of 5%, 15% relative risk reduction (RRR), and statistical test power of 80%. We performed TSA with a fixed-effects model when I 2 was less than 30%. Elsewise, random-effects model would be executed.

| Expression and survival analysis based on TCGA
The Tumor Immune Estimation Resource (TIMER) is an online database that incorporates expression profiles of 10,009 samples across 23 cancer types from TCGA (https://cistr ome.shiny apps.io/timer/). 22 We utilized TIMER to confirm the expression level of POU5F1 in various cancers. Gene Expression Profiling Interactive Analysis (GEPIA), another online analysis tool, contains survival and clinicopathological data extracted from various cancers based on TCGA (http:// gepia.cance r-pku.cn/). 23 Survival analyses of OS and DFS were executed by GEPIA to find the correlation between POU5F1 expression and the prognosis of various cancers.

| Specimens
Plasma specimens of 30 LIHC patients from Zhongnan Hospital of Wuhan Universitywere collected during July 2017 and October 2019 and stored at −80°C until use. LIHC patients were identified on the basis of their pathology reports. Meanwhile, 30 healthy people without hepatic diseases or abnormal liver biochemical outcomes were enrolled as controls. Our study was authorized by the Medical Ethics Committee of Zhongnan Hospital of Wuhan University.

PCR analysis
RNA was extracted from plasma using Total RNA Separate Extraction Kit (Bioteke, China) according to the manufacturer's instructions. NanoDrop 2000C was applied to assess the concentration and purity of RNA. ReverTra Ace qPCR RT Kit (Toyobo, Japan) was used to reverse transcript mRNA into complementary DNA (cDNA). Quantitative PCR (qPCR) was implemented using SYBR Green I UltraSYBR Mixture (CWBIO) on Bio-Rad CFX96 (Bio-Rad Laboratories). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the endogenous reference gene. The detailed sequences of each pair of primers are listed in Table S1. All experiments were repeated twice. The expression status of the target gene was assessed by the 2 −ΔCq method, in which ΔCq represents the value of the mean quantification cycle (Cq) of the target gene subtracting the mean Cq of the endogenous reference gene.

| Tumor infiltrating immune cell reckoning
CIBERSORT provides a deconvolution algorithm that is able to distinguish 22 kinds of tumor infiltrating immune cells (TIICs) from other cell types in tissues. 24 Expression profiles of 50 normal liver tissues and 374 LIHC tumor tissues were downloaded from TCGA database, and TIIC proportions of each sample were evaluated by R (version 3.6.2) on the basis of the CIBERSORT algorithm. Then, the TIIC proportions of normal liver tissues and LIHC tumor tissues were divided into two subgroups based on the median POU5F1 expression level and visualized through violin plots.

| Gene set enrichment analysis
Gene set enrichment analysis (GSEA) is a bioinformatics method that inspects the statistical significance of a priori defined sets of genes and verifies the differences between two biological states. 25 We divided TCGA LIHC samples into two phenotype subgroups on the grounds of the median expression level of POU5F1. Genes from the TCGA expression profiles were ranked in a list according to the degree of divergence between the high POU5F1 subgroup and the low POU5F1 subgroup through GSEA software 4.0. Then, Gene Ontology (GO) gene sets and Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets were analyzed to identify functional terms and pathways enriched in each phenotype subgroup. Gene set permutations were executed 1000 times for each analysis. The criteria of significantly enriched pathways were normalized p value < .05 and the absolute value of normalized enrichment score (NES) > 1.5.

| Enrichment analysis of POU5F1 coexpression genes
Co-expression genes of POU5F1 were screened out by R based on the expression profiles of TCGA. The cut-off line was set at p < .05, and the absolute value of Spearman correlation coefficient > 0.45. GO enrichment analysis and KEGG pathway analysis were performed using the R package "clus-terProfiler". p < .05 was taken as the statistically significant threshold.

| PPI network establishment and hub gene identification
We utilized the Search Tool for the Retrieval of Interacting Genes (STRING) database to establish a protein-protein interaction (PPI) network and discover the relationship among co-expression genes of POU5F1. The interaction score was set at 0.4 in the STRING database. Cytoscape was used to enhance the legibility of the PPI network on the basis of interaction data obtained from the STRING database. We considered genes that interacted directly with POU5F1 as central hub genes and those that directly interacted with the central hub genes as subordinate hub genes.

| Statistical analysis
All statistical analyses in this study were performed using Stata SE15 (Stata Corporation), SPSS 25.0 (SPSS Inc.), GraphPad Prism 8.0 (GraphPad Inc.), and R (version 3.6.2). Amalgamative HRs and relating 95% CIs in the meta-analysis were computed by Stata SE15. If the studies did not provide HRs or corresponding 95% CIs, these values were calculated by the equation in which P0 and P1 represented the survival rates of the decreased and elevated POU5F1 subgroups, respectively. The 95% CI was calculated through exp (lnHR ±1.96 × stderr), exp represented exponential, lnHR was natural logarithm of HR, and stderr meant standard error of HR. Several studies did not report the relevant data about survival rate in subgroups, and we utilized Engauge Digitizer Version 10.8 to collect representative data on Kaplan-Meier survival curves. Then, the extracted data were imported into a computation sheet obtained from Tierney et al. for estimation of HRs and 95% CIs. 26 Heterogeneity of enrolled cohorts was evaluated through Chi-square-based Q and I 2 analyses. We ran a metaanalysis with a fixed-effects model when the heterogeneity was acceptable (I 2 < 50% or p > .05). Otherwise, the randomeffects model was performed. Sensitivity analysis was conducted by sequentially expurgating each cohort to evaluate the stability of the amalgamative result. Potential publication biases were detected through Begg's and Egger's analyses.
Continuous variables with normal distribution are described as the mean ±standard deviation (SD). Median and inter-quartile ranges were used to describe abnormally distributed continuous variables. Student's t test and Mann-Whitney U tests were utilized for comparisons between two groups. Correlation analyses were conducted by Pearson and Spearman correlation tests. Chi-square test and Fisher's exact test were applied for evaluation of categorical variables. Cox proportional hazard regression was used for univariate and multivariate analyses. Odds ratios (ORs) were calculated by logistic regression. Receiver operating characteristic (ROC) curve was performed to assess the diagnostic values. The statistically significant threshold of the two-sided p value was set at .05.

| Literature search results and quality evaluation
The retrieval procedure is illustrated in Figure S1. The search of the databases obtained 1542 references. A total of 1205 articles remained after the exclusion of duplicates. About 909 records were excluded by scanning the titles and abstracts. 296 studies were evaluated by browsing the full-text. Eventually, 57 studies containing 7401 patients were enrolled in our study. [10][11][12][13][14][15][16][17][18]  are exhibited in Table S2. The quality assessments of the enrolled studies were implemented through the NOS, and 45 studies were rated as high-quality studies with comprehensive scores greater than 7 points (Table S3).

| Overall analysis of POU5F1 expression and cancer prognosis
Among the included 57 studies, a total of 5,485 subjects in 48 studies described the relationship between POU5F1 expression and OS, 1649 subjects in 14 studies for DFS, five studies with 1249 subjects for DSS and six studies involving 636 subjects for RFS. According to the meta-analyses, the heterogeneities were not distinct in these four kinds of prognosis analyses ( Figure 1A,C). Therefore, we capitalized on a fixed-effects model to calculate the amalgamative HRs and relating 95% CIs. The results showed that increased POU5F1 was correlated with inferior outcomes for OS (

| Subgroup analysis for OS and DFS
Subgroup analyses were implemented for OS and DFS to clarify the connection between POU5F1 expression and cancer type, analysis type, sample size, and detection method. Studies were defined as "other cancers" in the cancer type subgroup when there was only one enrolled study for each kind of cancer. As demonstrated in Figure 1B and Table 1, elevated expression of POU5F1 predicted poor prognosis of OS in bladder cancer, breast cancer, colorectal cancer, esophageal cancer, gastric cancer, LIHC, head and neck cancer, lung cancer and other cancers, including ovarian cancer, cervical cancer, neuroblastomas, and pancreatic cancer. However, the prognostic value of POU5F1 was not obvious in the overall survival of acute myeloid leukemia patients. Simultaneously, overexpression of POU5F1 was related to shorter DFS in head and neck cancer, breast cancer, LIHC, colorectal cancer, and other cancers, including lung cancer, gastric cancer, cervical cancer, and acute myeloid leukemia ( Figure 1D, Table S4). Furthermore, the subgroup category of analysis type, sample size, and detection method also indicated an observable relationship between a high level of POU5F1 and shorter OS and DFS.

| Correlation of POU5F1 and clinicopathological characteristics
To explore why elevated POU5F1 could lead to worse prognosis in various cancers, the correlations between POU5F1 status and neoplastic clinicopathological parameters were evaluated ( Table 2). The overexpression of POU5F1 was remarkably correlated with tumor size, TNM stage, tumor differentiation, tumor invasion depth, lymph node metastasis, distant metastasis, lymphovascular invasion, vascular invasion, tumor number, and tumor recurrence. Non-statistically significant results were found in age, gender, tumor encapsulation, liver cirrhosis, HBsAg, and smoking.

| Reliability of pooled prognostic results
TSA was implemented to assess the reliability of our meta-analysis results ( Figure S2). The heterogeneity of OS (I 2 = 4.40%), DFS (I 2 = 20.49%), and DSS (I 2 = 0.00%) was not obvious, so the fixed model was utilized to perform TSA. However, heterogeneity appeared in RFS (I 2 = 32.53%); thus, a random model was adopted. The accumulated Z-curve of OS crossed the traditional boundary, TSA boundary, and APIS, suggesting that the conclusion was significantly reliable. The cumulative Z-curves of DFS, DSS, and RFS crossed the conventional boundary and TSA boundary but did not reach the APIS, indicating that the current trials have obtained positive results, and more studies are required to support the results. Sensitivity analyses were executed to detect the stability of the conclusions about the prognostic value of POU5F1. No individual cohort could distinctly affect the pooled HRs of OS, DFS, DSS, or RFS, meaning the conclusions were credible ( Figure  S3). The underlying publication bias was appraised through Begg's and Egger's analyses. There was no potential publication bias found in OS, DFS, DSS, or RFS ( Figure S4).

| Expression and prognostic role of POU5F1 in various cancers
To further verify the expression level of POU5F1 in various cancers, TIMER was adopted to analyze the expression profiles from TCGA. As displayed in Figure 2A  survival analyses were carried out through GEPIA based on TCGA. Among the above-mentioned cancers, only LIHC showed statistically significant differences in both OS and DFS ( Figure 2B, Figure S5). Hence, LIHC was selected as the main target to explore the underlying functional role of POU5F1.

| Association between POU5F1 and clinicopathological variables of LIHC
The expression profiles and clinical characteristics of 374 LIHC patients were obtained from TCGA to probe the relationship between POU5F1 expression status and clinicopathological characteristics. As shown in Figure 3A-H, elevated POU5F1 was associated with tumor occurrence (p < .001), advanced histological grade (p = .016), stage (p = .025), tumor invasion depth (p = .019), and distant metastasis (p = .018). Logistic regression analysis indicated the expression of POU5F1 as a risk factor that was associated with poor prognostic clinicopathologic variables (

| Diagnostic value of POU5F1 in plasma
We detected the expression level of POU5F1 in plasma collected from 30 LIHC patients and 30 normal controls by qPCR to investigate the diagnostic value of POU5F1. The main clinical characteristics of the enrolled subjects are listed in Table 5. Significantly higher alanine aminotransferase (ALT) (p < .001), aspartate aminotransferase (AST) (p < .001), γ-glutamyl transferase (GGT) (p = .047), AFP (p < .001), and glucose (GLU) (p < .001) levels were observed in LIHC patients. In contrast, albumin (ALB) was much lower in LIHC patients than in normal controls (p < .001). The qPCR results revealed that POU5F1 was upregulated in the plasma of LIHC patients, which was consistent with the results from liver tissue samples based on TCGA ( Figure 3I). Moreover, elevated POU5F1 was associated with a high level of ALT in plasma (p < .001) ( Table 6). ROC analysis was utilized to assess the diagnostic value of POU5F1 in LIHC. As displayed in Figure 3J, the predictive validity of POU5F1  (AUC = 0.790, Se = 73.3%, Sp =80.0%, p < .001) was higher than that of AFP (AUC = 0.766, Se = 63.3%, Sp = 100.0%, p < .001). Encouragingly, the diagnostic validity was remarkably improved through the combination of POU5F1 and AFP (AUC = 0.902, Se = 83.3%, Sp = 80.0%, p < .001).

POU5F1 and TIICs
To inquire into the mechanism of POU5F1 involved in the pathological progression of LIHC, we analyzed the correlation between POU5F1 expression and 22 types of TIICs through the CIBERSORT algorithm on the basis of expression profiles from TCGA. As exhibited in Figure  4A, T cells CD4 memory resting (p = .019), T cells regulatory (p = .048), macrophage M1 (p = .012), and dendritic cells resting (p = .002) increased in the high POU5F1 group of normal liver tissues, while T cells follicular helper (p = .031) decreased. In LIHC tumor tissues, B cells memory (p = .001) and T cells follicular helper (p = 007) were enriched in the high POU5F1 group, and B cells naive (p < .001), monocytes (p = .004), and dendritic cells activated (p = .006) were increased in the low POU5F1 group ( Figure 4B). In addition, B cells memory (p < .001), T cells follicular helper (p = .039), and dendritic cells activated (p < .001) were positively related to POU5F1 in LIHC tumor tissues ( Figure 4C-E). The anomalous correlation between POU5F1 and dendritic cells activated might be partially explained by the limited data from dendritic cells activated. Negative correlations were observed in B cells naive (p < .001) and monocytes (p = .011) with POU5F1 ( Figure 4F, G).

| Identification of POU5F1related pathways
POU5F1-related signaling pathways were analyzed through GSEA to identify pathways that were differentially activated in LIHC between low and high POU5F1 expression phenotypes. GO terms enriched in the high POU5F1 phenotype mainly contained DNA replication, regulation of cell cycle G2M phase transition, signal transduction by p53 class mediator and so on. GO terms, including acute phase response and complement activation alternative pathway, were enriched in the low POU5F1 phenotype ( Figure 5A). Multiple cancer-related KEGG pathways were enriched in the high POU5F1 phenotype, such as bladder cancer, colorectal cancer, non-small cell lung cancer, and renal cell carcinoma. Several well-known cancer-related signaling pathways were also enriched in the high POU5F1 phenotype, including the MTOR signaling pathway, p53 signaling pathway, and WNT signaling pathway. The PPAR signaling pathway was enriched in the low POU5F1 phenotype ( Figure 5B).

| Enrichment analysis of POU5F1 coexpression genes
To explore genes that might potentially be associated with POU5F1, co-expression analysis was performed, and the expression status of the top 50 genes are displayed in Figure  5C. GO functional enrichment analysis indicated that these genes were enriched in cell proliferation-related terms, including chromosome segregation, nuclear division, organelle fission, and DNA replication ( Figure 5D). Cell cycle, spliceosome, p53 signaling pathway, pancreatic cancer, and DNA replication were the main signaling pathways in which these POU5F1 co-expression genes were enriched through KEGG pathway analysis ( Figure 5E).

| PPI network construction and hub gene recognition
A PPI network was constructed to reveal the intrinsic correlations among the POU5F1 co-expression genes. As exhibited in Figure 6A, the deeper color of each gene circle indicated an increased correlation coefficient with POU5F1. Analogously, a larger circle size indicated a smaller P value. Three genes (CBX3, CCHCR1, and NFYC) were found to be directly associated with POU5F1 and were defined as central hub genes. BARD1, ZNF692, IQCC, FBXL19, GPD2, and KAT2A had direct connections with the central hub genes and were regarded as subordinate hub genes for POU5F1. The expression status of the central hub genes and subordinate hub genes were all positively correlated with the expression level of POU5F1 in LIHC on the basis of TCGA ( Figure 6B). In addition to KAT2A, shorter OS of LIHC was found to be correlated with the overexpression of all the hub genes ( Figure 7A-I). Elevated expression of all the hub genes except NFYC indicated poor DFS of LIHC ( Figure 7J-R). In addition, the nine hub genes were all prominently upregulated in LIHC patients based on TCGA ( Figure 8A).

| DISCUSSION
POU5F1 has been studied for a long period of time as a wellknown CSC marker that participates in tumor invasion, differentiation, and recurrence. 67   on TCGA, which was consistent with the meta-analysis results. Interestingly, differences in both OS and DFS between the high POU5F1 group and the low POU5F1 group were observed only in LIHC on the basis of TCGA, indicating that POU5F1 played a unique and important role in the prognosis of LIHC. Furthermore, DNA replication, regulation of cell cycle G2M phase transition, bladder cancer, colorectal cancer, non-small cell lung cancer, renal cell carcinoma, MTOR signaling pathway, p53 signaling pathway, and WNT signaling pathway were the main GO and KEGG terms enriched in the high POU5F1 phenotype according to the GSEA. The GSEA results suggested that POU5F1 might participate in the pathological progression of LIHC and other cancers by promoting cell proliferation. Similar GO terms and KEGG pathways were found in the co-expression genes of POU5F1 and further validated the GSEA results.
Previous studies reported that TIICs could independently predict OS among cancer patients and reflect the status of lymph nodes. 75 Our study found that there was a prominent decrease in B cells naive and an increase in B cells memory in the high POU5F1 group of LIHC tumor tissues, hinting that the elevated POU5F1 might promote the transformation of B cells naive into B cells memory in LIHC. T cells follicular helper decreased in the high POU5F1 group in normal liver tissues, but increased in the high POU5F1 group in LIHC tumor tissues. The exact opposite results were found in dendritic cells activated. The typing and quantity conversion of TIICs in normal tissues and tumor tissues indicated the significant meanings of POU5F1 in regulating the tumor immune microenvironment of LIHC. The mechanism by which POU5F1 participates in the regulation of the tumor immune microenvironment still needs further study.
We identified upregulated POU5F1 as an independent prognostic factor for poor prognosis of LIHC through Cox regression, along with tumor stage, invasion depth, and distant metastasis. Overexpression of POU5F1 was related to high levels of ALT in plasma. Considering that a high concentration of ALT was indicative of liver cell destruction, we speculated that the overexpression of POU5F1 might be associated with hepatocellular necrosis or apoptosis. 76 In addition, although the diagnostic value of POU5F1 in LIHC was quite gratifying, the necessity of applying POU5F1 and AFP together in the diagnosis of LIHC to improve the diagnostic specificity needed to be emphasized, in view of POU5F1, was upregulated in a variety of cancers and might reduce the diagnostic specificity.
The molecular regulation mechanisms and pathways by which POU5F1 participates in LIHC have not been thoroughly studied. To further explore the role of POU5F1 in LIHC, we constructed a PPI network using co-expression genes of POU5F1 and identified hub genes that interacted with POU5F1, including CBX3, CCHCR1, NFYC, BARD1, ZNF692, IQCC, FBXL19, GPD2, and KAT2A. Based on related studies of hub genes, we visualized the pathways POU5F1 might play a role in LIHC ( Figure 8B). It has been reported that the transcription factor complex of POU5F1, SOX2, and KLF4 binds to the Nanog promoter to induce cellular reprogramming and cancer stemness. 77 EpICD translocates to the nucleus in a multiprotein complex and enhances the expression of POU5F1 by binding to the promoter of POU5F1. 78 CBX3 has been confirmed to promote cell cycle transition by inducing CDK1 and PCNA. 79 The elevated POU5F1 in LIHC may influence the expression of CBX3 and then activate the NF-Kβ and PI3K/Akt pathways through BARD1 and ZNF692. 80,81 The promotion effect of CCND2 on cell proliferation is regulated by NFYC and may also be affected by POU5F1. 82 The interaction between NFYC and KAT2A indicates that POU5F1 participates in tumor development through KAT2A-mediated histone H3 succinylation. 83 GPD2 promotes HuH-7 cell mitochondrial energy metabolism which may be regulated by NFYC and POU5F1. 84 As a central hub gene of POU5F1, CCHCR1 accelerates cell proliferation through EGFR. 85 FBXL19 induces Rac1 and Rac3 expression and inhibits apoptosis. 86 The relationship among FBXL19, CCHCR1, and POU5F1 needs further verification. In addition, loss of function of POU5F1 remarkably restrains propagation, metastasis, and aggression of cancer stem cells through inhibition of the PI3K/Akt pathway, from which we could expect POU5F1 to be an underlying target for cancer therapy. 87

| CONCLUSION
In summary, our study identified POU5F1 as a pan-cancer gene with significant prognostic value in various cancers, especially in LIHC. POU5F1 can serve as an independent prognostic factor for LIHC, and the combination of AFP and POU5F1 in plasma has prominent diagnostic validity for LIHC. POU5F1 may influence the progression of LIHC by regulating the tumor immune microenvironment and participating in cell proliferation-related pathways. Further research should be performed to verify the functional mechanism of POU5F1 in the pathogenesis of LIHC.