Literature search results and quality evaluation
The retrieving procedure was illustrated in Additional file 2. The databases search obtained 1542 references. 1205 articles were left after exclusion of duplicates. 909 records were excluded by scanning the titles and abstracts. 296 studies were evaluated through browsing the full-text. Eventually, 57 studies containing 7401 patients were enrolled into our study [10-18, 27-74]. 16 types of cancers were contained, including acute myeloid leukemia, bladder cancer, breast cancer, cervical cancer, colorectal cancer, esophageal cancer, gallbladder adenocarcinoma, gastric cancer, head and neck cancer, LIHC, lung cancer, neuroblastomas, ovarian cancer, pancreatic cancer, papillary renal cell carcinoma and prostate cancer. The expression levels of POU5F1 were detected by immunohistochemistry (IHC) in 44 studies, qPCR in 10 studies, and immunofluorescence (IF) in the remaining 3 studies. Essential features of enrolled studies were exhibited in Additional file 3. The quality assessments of enrolled studies were implemented through NOS and 45 studies were rated as high-quality studies with comprehensive scores greater than 7 points (Additional file 4).
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, 5 studies with 1249 subjects for DSS and 6 studies involved 636 subjects for RFS. According to the meta-analyses, the heterogeneities were not distinct in these 4 kinds of prognosis analyses (Fig. 1a, c). Therefore, we capitalized 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 (HR = 2.45, 95% CI = 2.22−2.71, P < 0.001), DFS (HR = 2.66, 95% CI = 2.22−3.19, P < 0.001), DSS (HR = 4.03, 95% CI = 2.70−6.01, P < 0.001) and RFS (HR = 2.59, 95% CI = 1.85−3.63, P < 0.001).
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, 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 Fig. 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. But the prognostic value of POU5F1 was not obvious in overall survival of acute myeloid leukemia. 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 (Fig. 1d, Additional file 5). Furthermore, the subgroup category of analysis type, sample size and detection method also indicated the observable relationship between 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 about age, gender, tumor encapsulation, liver cirrhosis, HBsAg and smoke.
Reliability of pooled prognostic results
TSA was implemented to assess the reliability of our meta-analysis results (Additional file 6). The heterogeneity of OS (I2 = 4.40%), DFS (I2 = 20.49%) and DSS (I2 = 0.00%) was not obvious, so the fixed model was utilized to perform TSA. Whereas, heterogeneity appeared in RFS (I2 = 32.53%), thus random model was adopted. The accumulated Z-curve of OS crossed traditional boundary, TSA boundary and APIS, suggesting the conclusion was significantly reliable. The cumulative Z-curve of DFS, DSS and RFS crossed the conventional boundary and TSA boundary but didn’t reach APIS, indicating the current trails have obtained positive results and more studies were 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 (Additional file 7). The underlying publication bias was appraised through Begg’s and Egger’s analysis. And there was no potential publication bias found in OS, DFS, DSS or RFS (Additional file 8).
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 Fig. 2a, POU5F1 was prominently upregulated in bladder urothelial carcinoma (BLCA); breast invasive carcinoma (BRCA); cholangiocarcinoma (CHOL); colon adenocarcinoma (COAD); head and neck squamous cell carcinoma (HNSC); kidney renal clear cell carcinoma (KIRC); kidney renal papillary cell carcinoma (KIRP); LIHC; rectum adenocarcinoma (READ); stomach adenocarcinoma (STAD); uterine corpus endometrial carcinoma (UCEC). Interestingly, downregulation of POU5F1 was only observed in kidney chromophobe (KICH). Further, survival analyses were carried out through GEPIA based on TCGA. Whereas, among the above-mentioned cancers, only LIHC showed statistically significant differences in both OS and DFS (Fig. 2b, Additional file 9). Hence, LIHC was selected as the main target to explore the underlying function role of POU5F1 played in.
Association between POU5F1 and clinicopathological variables of LIHC
The expression profiles and clinical characteristics of 374 LIHC patients were obtained from TCGA to probe into the relationship between POU5F1 expression status and clinicopathological characters. As shown in Fig. 3a-h, elevated POU5F1 was associated with tumor occurrence (P < 0.001), advanced histological grade (P = 0.016), stage (P = 0.025), tumor invasion depth (P = 0.019), and distant metastasis (P = 0.018). Logistic regression analysis indicated the expression of POU5F1 as a risk factor that was associated with poor prognostic clinicopathologic variables (Table 3). Increased POU5F1 was significantly correlated with tumor occurrence (OR = 65.63, P < 0.001), advanced stage (OR = 2.06, P = 0.007) and tumor invasion depth (OR = 2.00, P = 0.001). Besides, POU5F1 (HR = 1.64, P = 0.038) was identified as an independent risk factor for OS of LIHC through multivariate analysis, as were tumor stage (HR = 1.51, P < 0.001), tumor invasion depth (HR = 1.51, P < 0.001) and distant metastasis (HR = 3.73, P = 0.026) (Table 4).
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 characters of enrolled subjects were listed in Table 5. Significantly higher alanine aminotransferase (ALT) (P < 0.001), aspartate aminotransferase (AST) (P < 0.001), γ-glutamyl transferase (GGT) (P = 0.047), AFP (P < 0.001) and glucose (GLU) (P < 0.001) were observed in LIHC patients. On the contrary, albumin (ALB) was much lower in LIHC patients compared with normal controls (P < 0.001). The results of qPCR revealed POU5F1 was upregulated in plasma of LIHC patients, which was consistent with the results from liver tissue samples based on TCGA (Fig. 3i). Moreover, elevated POU5F1 was associated with high level of ALT in plasma (P < 0.001) (Table 6). ROC analysis was utilized to assess the diagnostic value of POU5F1 in LIHC. As displayed in Fig. 3j, the predictive validity of POU5F1 (AUC = 0.790, Se = 73.3%, Sp = 80.0%, P < 0.001) was higher than AFP (AUC = 0.766, Se = 63.3%, Sp = 100.0%, P < 0.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 < 0.001).
Relationship between POU5F1 and TIICs
In order to inquire into the mechanism of POU5F1 involved in the pathological progress of LIHC, we analyzed the correlation between POU5F1 expression and 22 types of TIICs through CIBERSORT algorithm on the basis of expression profiles from TCGA. As exhibited in Fig. 4a, T cells CD4 memory resting (P = 0.019), T cells regulatory (P = 0.048), macrophage M1 (P = 0.012) and dendritic cells resting (P = 0.002) increased in the high POU5F1 group of normal liver tissues, while T cells follicular helper (P = 0.031) decreased. In LIHC tumor tissues, B cells memory (P = 0.001) and T cells follicular helper (P = 0.007) were enriched in the high POU5F1 group, B cells naive (P < 0.001), monocytes (P = 0.004) and dendritic cells activated (P = 0.006) increased in the low POU5F1 group (Fig. 4b). Besides, B cells memory (P < 0.001), T cells follicular helper (P = 0.039) and dendritic cells activated (P < 0.001) were positively related with POU5F1 in LIHC tumor tissues (Fig. 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 < 0.001) and monocytes (P = 0.011) with POU5F1 (Fig. 4f, g).
Identification of POU5F1 related 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 high POU5F1 phenotype mainly contained DNA replication, regulation of cell cycle G2 M 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 low POU5F1 phenotype (Fig. 5a). Multiple cancers 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 cancers related signaling pathways were also enriched in the high POU5F1 phenotype, including MTOR signaling pathway, p53 signaling pathway and WNT signaling pathway. While, PPAR signaling pathway were enriched in the low POU5F1 phenotype (Fig. 5b).
Enrichment analysis of POU5F1 co-expression genes
To explore genes that might potentially associated with POU5F1, co-expression analysis was performed and the expression status of the top 50 genes were displayed in Fig. 5c. GO functional enrichment analysis indicated these genes were enriched in cell proliferation-related terms, including chromosome segregation, nuclear division, organelle fission, DNA replication and so on (Fig. 5d). Cell cycle, spliceosome, p53 signaling pathway, pancreatic cancer and DNA replication were the main signaling pathways in which these POU5F1 co-expression genes enriched through KEGG pathway analysis (Fig. 5e).
PPI network construction and hub genes recognition
PPI network was constructed to reveal the intrinsic correlations among the POU5F1 co-expression genes. As exhibited in Fig. 6a, deeper color of each gene circle indicated increased correlation coefficient with POU5F1. Analogously, larger circle size indicated smaller P value. 3 genes (CBX3, CCHCR1, 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 (Fig. 6b). In addition to KAT2A, shorter OS of LIHC was found to be correlated with the overexpression of all the hub genes. Elevated expression of all the hub genes except NFYC indicated poor DFS of LIHC (Fig. 7). Besides, the 9 hub genes were all prominently upregulated in LIHC patients based on TCGA (Fig. 8a).