{'fraction_non_alphanumeric': 0.03291982817455538, 'fraction_numerical': 0.015175237865831626, 'mean_word_length': 5.150617283950617, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of a novel histone phosphorylation prognostic signature in hepatocellular carcinoma based on bulk and single-cell RNA sequencing

Background Hepatocellular carcinoma (HCC) is the third leading cause of death in the world, characterized by high morbidity, poor prognosis and high mortality. Histone modifications regulate intracellular gene expression at the post-transcriptional level, and disturbances in the regulatory pattern of histone modifications at individual locus or across the genome can lead to tumorigenesis of HCC. In this study, we constructed a prognosis-related histone phosphorylation regulated (HPR) genes signature and elucidated whether HPR genes can predict overall survival in HCC patients. Methods Differentially expressed genes were screened using TCGA, ICGC and GEO databases, and a new risk signature was constructed by univariate Cox regression and Lasso regression analysis. Predictive nomograms were established by multivariate Cox regression of risk scores and clinical parameters, calibration curve and decision curve analysis were used to evaluate the models. The ssGSEA methods were used to determine the effect of risk scores on the tumor immune microenvironment. Data for HCC single-cell RNA sequencing (scRNA-seq) have been downloaded from Gene Expression Omnibus (GEO) to understand the role of HPR genes in tumorigenesis. Results Our analyses of nine HPR genes provided prognostic insights. Overall survival in the low-risk and high-risk groups was statistically higher, respectively (P<0.001). Cox regression analysis revealed that the risk score is a significant predictor of HCC outcomes (HR=2. 2.62, 95%CI: 1.248-5.514, P=0.011). In addition, a nomogram combining risk scores with TNM stages was constructed and tested from calibration curves and decision curves (AUC=0.780). MHC-class-I genes, iDCs, Macrophages, Tfh, Treg, Th2 were overexpressed in the high-risk group. Conclusion HPR genes risk score is closely related to the prognosis of HCC, tumor immune process and tumor cell progression.


Introduction
Hepatocellular carcinoma (HCC) is one of the most common aggressive malignant tumors. It ranks second in incidence and third in mortality among malignant tumors in the world, and is also the third leading cause of cancer-related death worldwide (1). HCC is a major concern facing the global today in health care, with nearly 800,000 new cases each year, its incidence is on the rise in worldwide (2)(3)(4). In a 2018 World Health Organization report, liver cancer was the second leading cause of cancer-related death after lung cancer (5,6). Despite recent breakthroughs in the treatment of viral hepatitis, the incidence of HCC is expected to continue to rise. With the rapid development of genomics, although in the same clinical stage, there are still great differences in the molecular characterization of tumor cells. The clinical efficacy and prognosis of treatment also vary widely according to the molecular typing of different types of patients. Various of the prognostic biomarkers commonly used in clinical practice still not satisfactory due to the varying levels of heterogeneity in HCC. HCC is mainly treated by surgery, combined with chemotherapy, radiotherapy and intervention therapy (7,8). Despite these multiple therapies, the prognosis of HCC patients is currently not favorable, and it is clinically necessary to actively construct a prognostic assessment model (9). The development of bioinformatics has greatly facilitated the screening of prognostic biomarkers in tumor patients. For example, Song et al. construct an endoplasmic reticulum stressrelated genes (ERGs) signature to predict the overall survival (OS) of patients with HCC (10). A prognostic model of lysine acetylation modifications (LARs) was developed to aid in individual prognostic monitoring and clinical decision-making in HCC (Sun et al.) (11). Fan et al. studied the expression profile and prognostic value of histone deacetylation in glioma for precise prediction of glioma prognosis (12). However, studies on the prognostic value of HPR in HCC are limited.
In recent years, histone modifications related to tumor formation and development have gradually attracted the attention of epigenetic researchers, and the increasing number of results indicate that tumors are often caused by histone modifications that affect their occurrence and prognosis (13). The epigenetic modifications offer ideal mechanisms for regulating gene expression and metabolic reprogramming. Histone phosphorylation is one of many ways of post-translational modification of histones. Histone phosphorylation is mainly involved in DNA damage repair, chromatin condensation, gene transcription and other cellular processes, and plays an important role in cell proliferation, differentiation, apoptosis and other processes. Aberrant histone phosphorylation has been detected in breast, prostate and colorectal cancer (14). In patients with chronic viral hepatitis, phosphatase overexpression is considered a key early risk factor for HCC. Aurora kinases are accountable for the phosphorylation of most serine/threonine residues in histone H3 (15). Aurora A is commonly overexpressed in HCC patients and is related with higher tumor grade. Activating NF-kB signaling through Aurora A reduces radiotherapy-induced apoptosis, thus preventing HCC from succumbing to radiation. A direct phosphorylation of histone H3 at Ser10 and Ser28 occurs after Aurora B interactions with histone H3, leading to chromosome instability and mitotic chromosome condensation (16,17). Aurora B is considered to be an independent molecular marker for predicting tumor growth and invasion in HCC (18). Several preclinical studies have shown that inhibition of Aurora B can inhibit tumor growth and induce apoptosis in different cancers (19,20). Current studies have shown that histone phosphorylation promotes the development of alcoholic liver disease. Ethanol and its metabolites increase histone H3 phosphorylation at serine 10 (H3S10) and serine 28 (H3S28) in primary rat hepatocytes by activating p38 mitogen-activated protein kinase (MAPK).
In this study, leveraging multi-omics data, we focused on identifying prognosis-related HPR genes and constructing an HPR-based genes signature, which was then validated for robustness with two external patient cohorts. In addition, the immune status of high-risk and low-risk groups was compared. We anticipate that this novel genetic signature may contribute to an independent prognostic indicator for HCC patients and facilitate the development of effective individualized therapy in the clinical setting.

Material and methods
Data screening and processing of bulk RNA-Seq data samples, along with details about clinical information. The International Cancer Genome Consortium (ICGC) data portal (https://dcc.icgc.org/) also provided clinical and mRNA expression data for the Japanese cohort of HCC (n = 231). At the same time, the dataset from the GEO website (https://www. ncbi.nlm.nih.gov/geo/) GSE14520 contains 221 United States cohort of HCC mRNA expression data and clinical information were also downloaded and used.
HPR-related genes were obtained from Gene Ontology Database (http://geneontology.org/) and Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/ msigdb/).After integration and deduplication, we established an optimized 47-gene HPR genes set. With the "limma" package in R, we identified HPR differentially expressed genes (DEGs) between HCC and normal tissue, based on FDR <0.05 and fold change >2. In parallel, univariate Cox regression was conducted to measure the correlation between each DEG and survival rate, and the genes showing a p value < 0.05 was considered a prognostic factor. Finally, differentially expressed HPGs associated with prognosis were obtained.

Acquisition and processing of single-Cell RNA-Seq data
The single-cell RNA-seq dataset GSE149614 for HCC was obtained from the GEO website, including 10 patients with primary tumors and 8 patients with normal liver tissue, and subsequent analyses were performed. Use the R package Seurat (21) (version:4.0.2) workflow to process single-cell data, normalize the data, find 2000 highly variable genes, perform PCA dimensionality reduction through highly variable genes, select the top 20 PCs and use the harmony package to multigroup data integration, and unified manifold approximation and projection (UMAP) is applied to perform dimensionality reduction and visual clustering of cells. By calculating mitochondrial gene content and the number of genes expressed in each cell, cells with mitochondrial gene expression less than 20% and the number of genes detected Between 500 and 7000 cells were retained, leaving 62737 cells after filtration. Finally, we used the FindAllMarkers function to screen for the fold change greater than 2(FC>2) and the expression ratio of the genes (Minpct = 0.5) to calculate the marker genes of the 25 subgroups, and then used the corrected p < 0.05 to screen the marker genes. The marker genes were matched with the cell groupings through the cell marker website, and the singleR package was used to assist in judging the cell groupings. Finally, 25 cell groups were identified.
We then extracted normal hepatocytes as well as cancer cells to construct normal cell-to-tumor cell trajectories using Monocle 2 (version: 2.10.1) (22). Results were visualized using the Plot Cell Trajectories and Plot Complex Cell Trajectories functions using highly variable genes for trajectory inference, and annotated with cell types and calculated cell states. Once the cellular change trajectories are defined, we examine the expression changes of candidate genes in the trajectories.
Construction and validation of the prognostic signature First, the differential genes were subjected to univariate cox regression to obtain HCC-related prognostic genes. In order to prevent overfitting, Lasso regression analysis was performed through the "glmnet" package in the R software to construct a prognostic signature. According to the corresponding genes in the Lasso Cox regression model, the coefficients and their expression levels yield the risk score for each patient: Score = e sum (expression of each gene × corresponding coefficient). A risk score is available for each patient, and based on the median, patients can be divided into two groups. Kaplan-Meier curve analysis assessed survival differences between the two groups. To visualize these two groups, we performed Principal Component Analysis (PCA) and t-distributed Stochastic Neighbor Embedding (t-SNE) using the "Rtsne" R package. Finally, the prognostic predictive value of the model was validated with the ICGC cohort and GEO cohort.

Visualization and assessment of prognostic value
We performed univariate and multivariate Cox regression analyses for multiple indicators that may affect survival outcomes, including risk score, sex, age, histological grade, AFP, fibrosis, and tumor stage. Based on multiple regression analysis, multiple predictors were integrated to construct a predictive nomogram model for OS. Calibration plots were used to compare the similarity between the 1-, 3-, and 5-year OS predicted using the nomogram with the idealized curves. The clinical utility of the predictive model was assessed using decision curve analysis.

Evaluation of immune status
To examine the link between risk scores and patient immune status, single-sample genomic enrichment analysis (ssGSEA) was performed to compare 16 distinct immune cell subsets and 13 immune function pathways in two risk cohorts using "GSVA" package in the R software. Boxplots created using the R packages "pRRophetic" and "ggplot2" visualize the results of the high-risk and low-risk cohorts.

Statistical analysis
Statistical analysis and data visualization were performed using R (version 4.1). Using the Wilcoxon test, differential genes were compared between HCC and matched non-tumor tissue. Kaplan-Meier analysis was used to estimate differences in survival between patients in different risk groups. Univariate and multivariate Cox regression analyses were used to screen for independent predictors of OS. For the "time-ROC" R package, the receiver operating characteristic (ROC) is used to evaluate the predicted probability of the risk model. In different risk groups, ssGSEA enrichment scores were compared between high and low risk groups by Mann-Whitney test. Pearson's correlation was used to calculate whether there was a significant correlation between risk scores and active immunization pathways. Statistical significance was determined by P-value and false discovery rate less than 0.05.

Result Patients' characteristics
In this study, the TCGA database contains 377 HCC patients, only 374 patients have mRNA expression data, and 370 of these 374 patients have complete survival time and survival status data. 231 HCC samples in the ICGC dataset and 221 HCC samples in the GEO dataset were included as the validation cohort and data collected in December 2021. The completed statistical characteristics of the cases are shown in Table 1, and the completed flow chart of this study is shown in Figure 1.  Schematic diagram of the overall design of the study.  (Figure 2A). First, we extracted the histone phosphorylation expression matrix, according to the screening criteria, a total of 23 genes showed differential expression in HCC compared with paracancerous tissues, including 22 upregulated and 1 down-regulated. The heat map ( Figure 2B) and volcano map ( Figure 2C) show the distribution of HPGs. Interactions between these prognostic genes are shown in correlation networks ( Figure 2D). Combining differential genes with patient survival information, univariate COX risk regression showed that 17 genes were associated with survival outcomes of HCC individuals ( Figure 2E).

Development and validation of histone phosphorylation gene signature
In order to ensure that an effective prognostic model can be predicted, the 17 prognostic genes of univariate analysis were compressed in the training group using LASSO regression to prevent overfitting, and 9 genes were obtained to build a model. of different groups of patients, and the low-risk group showed higher survival rates (p<0.001), a prognostic value exists for the risk score ( Figures 3B-D). ROC curves were widely used to assess the sensitivity and specificity of the model, and the accuracy of the histone phosphorylation gene signature in predicting 1-, 3-, and 5-year OS in our model was 0.742, 0.663, and 0.624 ( Figures 3E-G), respectively, indicates a good predictive power. In the TCGA cohort, ICGC cohort, and GEO cohort, all patients were divided into high and low risk groups using   Figures 4A-F). Dimensional reduction of patients by risk score was visualized, PCA and t-SNE analysis showed that the same risk group was clustered together (Figures 4G-L).

Construction of the nomogram
To further confirm whether risk score can serve as an independent prognostic factor for HCC, we performed multivariate and univariate Cox regression analysis together with possible clinical indicators, including TNM stage, gender, age, histological grade, AFP, fibrosis and risk score. As shown in Figure 5A-F, based on multivariate Cox regression analysis, risk score and clinical stage emerged as independent prognostic factors in the TGCA cohort (TCGA: HR = 2.62, 95% CI = 1.25-5.51, p =0.011). In order to obtain individualized assessment for each patient, we constructed nomograms of the variables obtained from the multivariate prognostic analysis, and display the constructed nomogram on the calibration curve and the ideal model can achieve similar results (Figures 5G-H). By combining the two independent prognostic factors, risk score and clinical stage, a significant improvement in mortality at 1-, 3-, or 5-year OS was seen (AUC=0.78) ( Figure 5I). Decision Curve Analysis (DCA) has also shown that predictive models have good clinical utility ( Figure 5J).

Correlation between risk score and clinical characteristics
We further investigated the clinical utility of the HPR risk score, examining the relationship between the risk score and clinicopathological features (including tissue grade, TNM stage, gender, and age). Figures 6A-G depict the distribution of risk scores between groups by clinical characteristics. Based on these data, we found that risk scores were independent of age and sex (p>0.05) (Figures 6A, B, E, F), whereas higher risk scores were found in patients with advanced TNM stage or higher pathological grade, and the ICGC cohort came to the same conclusion. In addition, Figures 6H-K show differences in survival between high and low risk groups at the same stage and grade. Among patients of the same grade and stage, survival analysis showed that the high-risk groups all had lower survival times. It indicated that the progression of HCC was related to HPR.

Immune status between two risk groups
In different risk groups, we measured enrichment scores for 16 immune cells and 13 immune functions using the ssGSEA algorithm. According to the results, antigen-presenting cells such as aDC, Th2, Treg, and macrophages were enriched in the high-risk group (Figures 7A, B). In addition, immune function levels including MHC class I were also enriched in the high-risk group, while NK cells, B cells, Type-II-IFN response, and Type-I-IFN response were enriched in the lowrisk group (Figures 7C-E). Furthermore, the presence of immune cells in the tumor microenvironment was positively associated with the risk score (including Tfh cells, aDC cells, Treg cells, Th2 cells and macrophages) and with immune function (MHC_class_I). However, NK cells, B cells, Type-II-IFN response, and Type-I-IFN was negatively correlated with the risk score shown in the lollipop plot ( Figure 7F). This difference in immune infiltration may account for poor prognosis in highrisk patients.

Validation of differential expression of prognostic genes in single-cell RNA-Seq data
After quality control and normalization of the dataset GSE136103, data for 21,560 genes and 62737 cells were retained. Using the umap algorithm, 62737 cells were divided into 25 clusters. We obtained human cell marker genes from CellMarker (http://biocc.hrbmu.edu.cn/CellMarker/), and based on the marker genes and the R package "singleR" to help mark cell grouping, we annotated these 25 clusters  (Figures 8A-E). We then examined the expression of 9 model genes in single cells and found that CDK1, AURKB, CCNA2, and CHEK1 are mainly expressed in CD4 cells, while PRKAA2, CDK5, and PRKCD are mainly expressed in tumor cells ( Figure 9A). To identify the developmental relationship between hepatocytes and malignant cells, we executed single-cell trajectory analysis of scRNA-seq data using Monocle. We therefore defined the hepatocyte cluster as the root and the tumor cell cluster as the end of the trajectory (Figures 9B-D). As cells transition from one state to another, some genes are silenced while others become newly active. We found that PRKCD, PRKAA2, CDK5, and CDK1 genes were increased in the progression from hepatocytes to tumor cells ( Figure 9E).

Discussion
HCC is a common malignancy, and the widely used TNM staging system is often ineffective in predicting prognosis. Biomarkers that predict patients' survival rates are urgently needed for HCC (23,24). Various types of histone phosphorylation are involved in gene transcription, DNA repair, apoptosis and chromosome condensation, and play a crucial role in metabolic regulation and tumorigenesis (25). Here, using multiomics HCC data from multiple cohorts, we performed an extensive bioinformatics analysis to build a risk model consisting of 9 HPR genes that were significantly associated with prognosis in HCC patients. By exploring the relationship between HPR score and clinical parameters in HCC patients, it was demonstrated that HPR was associated with high stage and grade, and immune infiltration was also analyzed to further investigate the clinical significance of prognosis. Finally, the expression of HPR genes in single-cell sequencing was analyzed, and the results also confirmed that most of HPR genes were highly expressed in tumor cells. This HPR risk model not only provides more intuitive information for clinicians to personalize the treatment of HCC in the future, but also provides a molecular basis for HPR. In this study, a prognostic risk model for 9 genes INCENP, NEK11, AURKB, CCNA2, PRKAA2, CHEK1, CDK5, PRKCD,   (26)(27)(28). Aurora-B is a chromosomal passenger protein that forms a complex with INCENP and survivin to regulate stable bipolar spindle-kinetochore attachment during mitosis, chromosome segregation and cytokinesis. Demonstrating that AURKB has been shown to be closely associated with liver cancer, its upregulation plays a key role in promoting hyper polyploidization, and an increase in AURKB phosphorylation was detected on intermediates during cytokinesis, leading to hyperpolyploidization (29-31). They are serine/threonine kinases that are collectively known as the NIMA-related kinase (NEK) family. Members of the family include NEK1-NEK11. NEK11 is integral to the cell cycle and microtubule formation (32). The expression of PRKAA2 helps activate autophagy and chemoresistance. PRKAA2 plays a key role in regulating autophagy and 5-FU resistance in gastric cancer (33). CHEK1 is responsible for the control of G1/S, S and G2/M checkpoints, and PRKCD inhibits autophagy. phagocytosis to promote the proliferation, invasion  and metastasis of hepatoma cells (34,35). The cyclin family gene CCNA2 has been shown to be tumor-promoting in multiple solid tumor types, and its expression is differential between many types of cancer and normal tissues. Cancer of multiple types may be affected by CCNA2 (36-38). The immune microenvironment plays an important role in tumorigenesis and has been widely recognized. The results of ssGSEA analysis showed that the high-risk group had higher infiltration levels of Treg, iDCs, macrophages and Th2 cells, and higher expression levels of MHC-class-I genes. Fu et al. showed that increased macrophages and Treg cells in HCC patients are associated with poor prognosis, and this is consistent with our findings in the high-risk group (39). Shankaran et al. elucidated that IFNg and lymphocytes in the immune system can cooperate to suppress tumor development, consistent with our results with elevated B cells, Type-II-IFN response in the low-risk group (40). The results demonstrate the relationship between HPR and immunity, emphaizing the critical role of immunotherapy in high-risk HCC patients.
scRNA-seq has become an important tool for studying tumor heterogeneity and tumor cell evolution in a variety of cancers (41). Here, we used HCC scRNA-seq data from the GEO database to explore trajectory analysis from normal liver cells to tumor cells. INCENP and NEK11 genes were expressed at a low level in the single-cell transcriptome due to the limitation of single-cell technology. CDK1, AURKB, CCNA2, and CHEK1 were highly expressed in immune cells, suggesting that histone phosphorylation is closely related to the tumor immune environment. Cell trajectories showed that the expression of PRKCD, PRKAA2, CDK5, and CDK1 genes increased from hepatocytes to tumor cells. These findings identify the nine genes that play important roles in immune cell infiltration and in the development of normal cells into tumor cells, providing a basis for subsequent molecular studies.
This study may have several unavoidable limitations. First, we constructed a differentially expressed HPR-based prognostic model to predict survival in HCC patients based on data from the TCGA database. Although two external datasets were used as validation, it still needs to be validated in real clinical samples. Second, differentially expressed HPR genes need to be validated in molecular experiments. The last major limitation is the lack of in-depth mechanism research. Bioinformatics provides an analytical basis for screening genes, but how these genes cause histone phosphorylation needs to be further improved by in vitro/in vivo experiments.

Conclusion
HCC is a common and high-burden tumor. We elaborated the r ole of H PGs in t he p rognosis and immu ne microenvironment of HCC. HP risk score has been identified as an important prognostic factor for HCC, improving the predictive power of the TNM staging system. Collectively, these findings provide new insights into the prognostic assessment and treatment of HCC.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.