A cohort of 453 patients from The Cancer Genome Atlas (TCGA)-LUAD was initially included for analysis. Each patient possessed comprehensive gene expression profiles along with corresponding clinical information. All data were obtained using the R package “TCGAbiolinks” (Colaprico et al., 2016). In addition, several cohorts released on Gene Expression Omnibus (GEO) database with overall survival (OS) time were also enrolled, including GSE31210, GSE68465, GSE72094, GSE11117, GSE11969, GSE13213, GSE = 42,127, GSE19188, GSE63459, GSE29016. GSE40791 and GSE31547 contains the gene expression matrix of both normal and tumor sample, and also be employed to compare the different expression. Importantly, patients with an overall survival (OS) time of less than 1 month were excluded to mitigate potential bias. For the TCGA-LUAD gene expression profile, genes with zero expression in more than 10% of samples were also excluded. The count data were converted to transcripts per kilobase million (TPM) values, followed by a log2(TPM+1) transformation for subsequent analysis (Lu et al., 2019). All the gene expression data ranges from 0 to 20 after scale. All the GEO cohorts can be downloaded from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Univariate Cox analysis was employed to identify prognostic genes in the TCGA-LUAD and three GEO cohorts, selecting those with a hazard ratio (HR) greater than 1.5 and a p-value less than 0.05. Subsequently, a Venn diagram was utilized to display the consensus prognostic genes for further analysis. Kaplan-Meier curves were generated to compare overall survival (OS) using the log-rank test. A meta-analysis based on the HR and 95% confidence intervals (95% CI) was calculated to uncover the overall prognostic value.

We calculated the correlations among genes with selected gene by Pearson correlation test, and the signaling enrichment of the top 200 genes were performed by MetaScape (http://metascape.org/) (Zhou et al., 2019). We evaluated the activated pathways using fast gene set enrichment analysis (fgsea, https://github.com/ctlab/FGSEA/). Initially, GSEA was conducted by ranking the input molecular readouts, followed by calculating the pathway enrichment score through a running-sum statistic. This method increases the score if a feature falls into the target pathway and decreases it otherwise. The final score represents the maximum deviation from zero observed during the random walk, normalized by computing the z-score of the estimate compared to a null distribution derived from random permutations. CancerSEA (Yuan et al., 2019) (http://biocc.hrbmu.edu.cn/CancerSEA) is the first dedicated database designed to comprehensively resolve the distinct functional states of cancer cells at the single-cell level. It provides a cancer single-cell functional state atlas encompassing 14 functional states—stemness, invasion, metastasis, proliferation, EMT, angiogenesis, apoptosis, cell cycle, differentiation, DNA damage, DNA repair, hypoxia, inflammation, and quiescence—across 41,900 cancer single cells from 25 cancer types. We downloaded these gene sets and utilized the z-score algorithm in the R package GSVA to calculate the functional status of the 14 gene sets, with the values for each gene set expressed as z-scores. Pearson correlation analysis was then employed to determine the statistical correlation of genes with each gene-based z-score.

We compared the different gene expression among wild type and mutated samples in TCGA-LUAD cohort by TIMER 2.0 (http://timer.cistrome.org/). The correlation between expression of selected gene and other mutated genes were also calculated by Pearson correlation analysis, and further validated in GSE26939 and GSE72094. Correlation between selected gene and immunocyte infiltration was also evaluated by TIMER 2.0, we also assessed the distribution of high and low selected gene expression TCGA-LUAD sample among six immune subtypes, including wound healing (C1), IFN-γ dominant (C2), inflammatory (C3), lymphocyte depleted (C4), immunologically quiet (C5), and TGF-β dominant (C6) (Thorsson et al., 2018). The potential response of immunotherapy and chemotherapy evaluated by BEST (https://rookieutopia.com/) (Zaoqu et al., 2023). IMvigor210 cohort contains 348 patients that received the therapy of PD-L1 blockade with atezolizumab in metastatic urothelial cancer (Mariathasan et al., 2018), the clinical information and gene expression data can be accessed from http://research-pub.gene.com/IMvigor210CoreBiologies/. The data of Wolf 2021 cohort can be assessed from GSE173839.

To observe the protein level and localization of specific gene, we checked the Human Protein Atlas (HPA, https://www.proteinatlas.org/) website, of which provide the picture of immunofluorescent staining. To further explore the distribution of specific gene in LUAD, we also employed the data from spatial transcriptomics and single-cell sequencing. We acquired the data about samples from patients with brain metastasis of non-small cell lung carcinoma (GSE179572) (Sudmeier et al., 2022), to evaluate the spatial distribution of specific gene. To accurately assess the cellular composition of each spot on the 10x Visium slides, we employed deconvolution analysis. This method leverages spatial transcriptomics and single-cell transcriptomics data, with particular consideration given to the specific cancer type. Based on the preceding deconvolution results, we calculated the predominant cell type in each microregion and visualized the highest cellular content in each microregion using the SpatialDimPlot function from the Seurat package (Supplementary Table S1). The SpatialFeaturePlot function from the Seurat package was utilized to visualize the gene expression landscape in each microregion. Spearman correlation analysis was performed to calculate the correlations between cellular contents across all spots and between cellular content and gene expression. The results were visualized using the linkET package. The data of single-cell profiling of advanced non-small cell lung cancer (GSE148071) (Wu et al., 2021) was also download for the further analysis. The analysis and virilization of single-cell data based on TISCH2 (http://tisch.comp-genomics.org/) (Han et al., 2023). The cell type of malignant and alveolar was collected from the source data, while the marker genes of immunocytes provided in Supplementary Table S2.

All the statistical analyses were performed by the R version 4.2.2. Student’s t-test was applied to compare two groups if the data is normally distributed for continuous data, otherwise Wilcoxon rank-sum test will be used. For the comparison of continuous data among more than two groups, Kruskal–Wallis test was utilized. As for categorical data, Chi-square test and Fisher’s exact test was conducted. Pearson correlation coefficient analysis was used to calculate the correlation of HJURP with other types of data. The receiver operating characteristic (ROC) area under the curve (AUC) was performed to access the stability of prediction. To find out the independent risk factors, univariate analysis and multivariate analysis were performed. p < 0.05 was considered statistically significant.

To identify prognostic genes, we employed univariate Cox regression analysis in four cohorts: TCGA-LUAD, GSE31210, GSE68465, and GSE72094. We assessed the prognostic predictive power of all genes, applying a predefined threshold (HR > 1.5, p < 0.01). This analysis revealed 474 risk genes in the LUAD cohort (Figure 1A), 968 risk genes in the GSE31210 cohort (Figure 1B), 121 risk genes in the GSE68465 cohort (Figure 1C), and 722 risk genes in the GSE72094 cohort (Figure 1D). Subsequently, we performed a merged analysis of the risk genes identified across the four cohorts by Venn plot (Figure 1E), uncovering 19 genes consistently associated with accelerated tumor progression. For all the 19 genes, HJURP shows the secondary rank with its expression fold change of 3.71 compared between tumor and normal samples, only less than KIF14, and for the expression level, HJURP have the higher level than KIF14 (Supplementary Figure S1). Holliday junction recognition protein (HJURP) is a key molecular chaperone for centromere protein A (CENP-A), which is essential for chromosome separation during mitosis and cell cycle regulation. Subsequently pan-cancer analysis, we also revealed that high HJURP expression in 18 types of tumors, including LUAD (all p < 0.05, Figure 1F).

Several studies based on clinical tumor samples also confirmed the risk of HJURP to the tumorigenesis or development of cholangiocarcinoma (Yang et al., 2022), colorectal cancer (Kang et al., 2020), pancreatic cancer (Wang et al., 2020), triple-negative breast cancer (Mao et al., 2022) and hepatocellular carcinoma (Chen et al., 2018). While HJURP has indeed been reported in other types of cancers, its role in LUAD remains underexplored, which presents a significant opportunity for novel discoveries. Our decision to focus on HJURP stems from its well-documented function in maintaining chromosomal stability through the centromeric loading of CENP-A, a key process that is often dysregulated in various cancers, including lung cancer.

We utilized Kaplan-Meier survival analysis to illustrate the predictive power of the HJURP gene across different cohorts. In the TCGA-LUAD cohort, patients with high HJURP expression had a 1.93-fold higher risk of overall survival (OS) death compared to those with low HJURP expression (95% CI: 1.321–2.815, p < 0.001, Figure 2A). In the GSE31210 cohort, the high HJURP expression group had a 2.75-fold higher risk of OS (95% CI: 1.319–5.735, p = 0.007, Figure 2B). In the GSE68654 cohort, high HJURP expression was associated with a 1.57-fold higher risk of OS (95% CI: 1.215–2.038, p < 0.001, Figure 2C). In the GSE72094 cohort, high HJURP expression corresponded to a 2.2-fold higher risk of OS (95% CI: 1.485–3.27, p < 0.001, Figure 2D).

Furthermore, we analyzed the impact of HJURP gene expression on various clinical outcomes within the TCGA-LUAD cohort. High HJURP expression was indicative of poorer disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI). Through a meta-analysis, we integrated the prognostic predictive power of HJURP across different cohorts and found that patients with high HJURP expression had a 1.87-fold higher risk of adverse prognostic outcomes compared to those with low expression (Figure 2E). These findings collectively demonstrate that high HJURP expression is significantly associated with decreased survival rates in multiple cohorts.

In the TCGA cohort, HJURP expression was significantly elevated in tumors, as compared between paired tumor and adjacent normal samples (p < 0.001, Figure 3A). This finding was further validated in the GSE40791 and GSE131547 cohorts (p < 0.001, Figure 3B). In clinical phenotype-related analyses, we found that HJURP expression was significantly higher in male patients compared to female patients (TCGA-LUAD: p = 0.002, GSE31210: p < 0.001, GSE68654: p = 0.012, Figure 3C). Moreover, HJURP expression was closely related to smoking status, with higher expression observed in current and ever smokers with lung cancer (TCGA-LUAD: p < 0.001, GSE31210: p < 0.001, Figure 3D), suggesting that smoking may influence tumorigenesis by upregulating HJURP expression. Additionally, we found that Asians had higher HJURP expression compared to other ethnicities (GSE68654: p < 0.001, Figure 3E).

HJURP expression increased with advanced tumor stage (TCGA-LUAD: p < 0.001, GSE31210: p < 0.001, GSE68654: p = 0.012, GSE72094: p = 0.012, Figure 4A). Compared to well-differentiated LUAD, poorly differentiated samples exhibited significantly higher HJURP expression (p < 0.001, Figure 4B). Additionally, we observed that patients with progressive disease in the TCGA-LUAD cohort had the highest HJURP expression in their tumor tissues (p = 0.036, Figure 4C). These results indicate that HJURP is closely associated with tumor progression, with males, smokers, and Asians being high-risk groups for HJURP-influenced LUAD development.

In further research, we conducted multivariate Cox regression analyses in different cohorts to eliminate confounding factors affecting HJURP’s prediction of LUAD prognosis (Table 1). In the TCGA-LUAD cohort, age over 70 years (p < 0.01), tumor stage (all other stages vs Stage I: p < 0.001), and HJURP expression (p < 0.001) were independent prognostic factors. In the GSE31210 cohort, age (p = 0.014), tumor stage (Stage II vs Stage I: p < 0.001), and HJURP expression (p = 0.042) were independent prognostic factors. In the GSE68465 cohort, age (p = 0.004), tumor grade (well vs moderate differentiation, p = 0.029), tumor T stage (T3 vs T1: p = 0.003, T4 vs T1: p = 0.004), tumor N stage (N1 vs N0: p < 0.001, N2 vs N0: p < 0.001), and HJURP expression (p < 0.001) were independent prognostic factors. In the GSE72094 cohort, gender (p = 0.003), tumor stage (Stage II vs Stage I: p = 0.004, Stage III vs Stage I: p < 0.001, Stage IV vs Stage I: p = 0.013), and HJURP expression (p < 0.001) were independent prognostic factors. Overall, multivariate regression analysis further confirmed HJURP expression as a predictor of LUAD prognosis, its high expression being a risk factor for poor prognosis.

With the immunofluorescent staining picture provided by HPA, we observed that HJURP mostly localized in the nucleoplasm and nucleoli (Figure 4D), which can reflect its potential function in the regulation of genetic alteration. Using spatial transcriptomics data, we observed that in LUAD tumor tissues (Figure 5A), the regions of HJURP expression (Figure 5B) correspond to tumor areas rather than regions populated by immune cells (Figure 5C). This indicates that HJURP is predominantly expressed in tumor cells. Additionally, we found that HJURP expression levels positively correlate with tumor cell density and negatively correlate with plasma cells, macrophages, endothelial cells, and fibroblasts (Figure 5D). In single-cell sequencing data based on lung cancer tissues, HJURP was significantly more highly expressed in malignant cells compared to immune cells and stromal cells (p < 0.001, Figure 5E). Specifically, we observed that regions with high HJURP expression in single-cell sequencing data highly overlap with regions of malignant cells (Figures 5F,G). Among cells with positive HJURP expression, malignant cells accounted for 80.6% (Figure 5H).

We calculated the expression correlation between HJURP and over 20,000 other genes (Figure 6A). We selected the top 200 genes with the highest expression correlation for biological function enrichment analysis. We found that HJURP may influence DNA metabolic processes, mitotic cell cycle processes, cell cycle phase transitions, and the cell cycle (Figure 6B). Using the FGSEA algorithm, we assessed the significantly different activation levels of signaling pathways between patients with high and low HJURP expression. In tumors of patients with high HJURP expression, cell cycle-related pathways were significantly activated, including MYC targets, G2M checkpoint, mitotic spindle, and E2F targets, as well as DNA repair, unfolded protein response, and glycolysis (Figure 6C). In further studies, we calculated the activation levels of 14 tumor development-related signaling pathways and evaluated their correlation with HJURP expression (Figure 6D). We found that HJURP was significantly positively correlated with the activation of cell cycle (R = 0.93, p < 0.001), DNA damage (R = 0.73, p < 0.001), DNA repair (R = 0.73, p < 0.001), cell proliferation (R = 0.44, p < 0.001), hypoxia (R = 0.26, p < 0.001), invasion (R = 0.28, p < 0.001), and metastasis (R = 0.17, p < 0.001) pathways, and showed a negative correlation with tumor differentiation (R = −0.16, p < 0.001) and stemness (R = −0.22, p < 0.001). In summary, HJURP may promote tumor progression in LUAD by influencing pathways related to the cell cycle and DNA damage repair.

Tumor development and progression are closely related to gene mutations and immune cell infiltration. Mutation of HJURP gene can alter its expression, with mutated HJURP exhibiting higher expression levels (p = 0.094, Figure 7A). HJURP mutation are also associated with poorer prognosis (p = 0.042, Figure 7C), likely due to its increased expression levels. Additionally, we found that HJURP expression is closely correlated with a series of gene mutations. Patients with high HJURP expression showed higher mutation frequencies in TP53, CSMD3, TTN, and MUC16 genes, while the mutation frequency of the EGFR gene was lower (p < 0.001, Figure 7B). These findings were validated in external cohorts, where HJURP expression was significantly elevated in TP53-mutated samples (all p < 0.001, Figure 7D) and decreased in EGFR-mutated samples (p < 0.001 for GSE72094, p = 0.55 for GSE26939, Figure 7E).

Regarding the correlation between HJURP expression and immune cell infiltration, our analysis revealed that high HJURP expression was significantly positively correlated with the infiltration of CD4-positive T cells, such as Th2 and Th1 cells, as well as M1 macrophages, while it was negatively correlated with M2 macrophages (Figure 7F). These results suggest that high HJURP expression may be associated with a pro-inflammatory state in the tumor microenvironment. Further analysis of the distribution of high and low HJURP expression groups in six immune subtypes showed that LUAD patients with high HJURP expression were more likely to belong to the IFN-γ dominant (C2) subtype (Figure 7G), which exhibited a high proliferation rate that may override an evolving type I immune response (Thorsson et al., 2018). HJURP expression also showed a positive correlation with IFN-γ (R = 0.33, p < 0.001, Supplementary Figure S2).

Based on previous research findings, we observed that HJURP expression is positively correlated with the activation of cell proliferation, DNA damage repair pathways, and certain immune cell infiltration and immune response activation. Therefore, we further analyzed potential precision therapy strategies influenced by HJURP expression. We examined the correlation between a series of immune checkpoint genes and HJURP expression and found a significant positive correlation between HJURP expression and immune checkpoints such as PD-L1, PD-L2, IDO1, and MICB across multiple cohorts, while the correlation with PD1 expression was less pronounced (Figure 8A). Further analysis revealed that in patients responding to anti-PD-L1 treatment, HJURP expression was significantly higher compared to non-responders (IMvigor210 cohort 2018: p < 0.001, Figure 8B; Wolf cohort 2021: p = 0.0032; Figure 8C), which was not the case for anti-PD1 treatment (Supplementary Figure S3).

Regarding potential chemotherapy treatments, we identified from the GDSC1 database that commonly used drugs such as mitomycin, and temozolomide may be effective against tumors with high HJURP expression (Figure 8D). Additionally, from the CTRP database, we found that sorafenib could potentially have therapeutic effects (Figure 8E). These three common anti-tumor drugs already be applicated in the clinical treatment of LUAD, while others predicted agents should be further validated.