This research analyzed three LIHC specimens (GSM3064824, GSM3064820, and GSM3064823) sourced from the GSE112271 dataset at the resolution of single cells. Furthermore, we merged RNA sequencing data alongside clinical information derived from the TCGA-LIHC dataset. In order to create and validate the diagnostic model, several datasets were employed, comprising TCGA-LIHC, GSE112790, and GSE102451. The clinical prognostic data of 240 primary liver cancer samples from Japan in the IGCG database were also included in this study.

The NMF algorithm was utilized to derive coefficients of biological significance from the gene expression matrix, arranging genes and samples to emphasize the structural properties of the data and aid in classification (19). Differential expression analysis for clusters A and B was carried out utilizing the ‘Limma’ R package, with the parameters established at |logFC| > 0.5 and an adjusted p-value lower than 0.05. Following this, the ‘NMF’ R package was employed to group all samples according to the differentially expressed genes (DEGs) found in the subclusters, with the goal of uncovering potential molecular subtypes. The ‘brunet’ algorithm was performed for 100 iterations for each specified value, varying from 2 to 10 clusters. The ideal number of clusters was identified by assessing cophenetic correlation, dispersion, and silhouette width (20). Furthermore, the Limma package in R (version 3.40.2) was used to analyze mRNA differential expression between malignant and adjacent non-malignant tissues in the TCGA-PRAD dataset.

To validate the reliability of the immune score results, we utilized the immunedeconv R package (21). An extensive assessment of every algorithm was conducted, highlighting their distinct benefits. The XCELL approach was chosen due to its ability to evaluate a wider variety of immune cell categories (22).

The training phase employed the TCGA-LIHC dataset, while validation was conducted using the GSE112790 and GSE102451 datasets. Model combinations were evaluated based on their area under the curve (AUC) values, identifying the optimal model as the one with the highest average AUC. Receiver operating characteristic (ROC) curve analysis was performed using the pROC package [1.18.0], with results visualized via ggplot2 [3.3.6].

The Limma package (version 3.40.2) in R was utilized to explore the differential expression of mRNA within the TCGA-LIHC dataset. We set thresholds for identifying differentially expressed mRNAs at “P < 0.05 and log2 (fold change) > 1.3 or log2 (fold change) < -1.3.” To further investigate the roles of target genes in carcinogenesis, we employed the ClusterProfiler package to analyze Gene Ontology (GO) functions and to enrich Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Relevant pathway-associated genes were compiled and analyzed using the GSVA package in R, and single-sample gene set enrichment analysis (ssGSEA) was conducted with the method parameter set to ‘ssgsea’. Additionally, we evaluated the correlation between gene expression and pathway scores using Spearman correlation analysis.

After initially baking the paraffin slices, they were immersed in two xylene baths, each lasting 15 minutes. This was followed by a series of consecutive soaks in absolute ethanol, 95% ethanol, 85% ethanol, 75% ethanol, and finally distilled water, with each solution applied for 5 minutes. The sections were then transferred to a retrieval box containing an alkaline antigen retrieval solution (pH 9.0 EDTA) and heated in a pressure cooker for 2 minutes. After allowing the sections to cool naturally, they were washed three times with PBS (pH 7.4) for 5 minutes each, with gentle stirring. Subsequently, the sections were treated with a 3% hydrogen peroxide solution for 15 minutes at room temperature in the dark to inhibit endogenous peroxidase activity. To ensure uniform tissue coverage, a blocking solution was applied, and the sections were blocked for 30 minutes at room temperature. Following this, diluted S100A10 antibody (11250-1-AP) and POU5F1 antibody (60242-1-Ig) were added and incubated overnight at 4°C. The next day, the sections were washed three times with PBS for 5 minutes each. After gently shaking to remove excess liquid, a poly-HRP secondary antibody corresponding to the primary antibody species was added dropwise and incubated in the dark at room temperature for 10 to 20 minutes. The final score of the staining result was calculated by multiplying the staining intensity by the staining range. The staining range was categorized as 0%-25% for 1 point, 26%-50% for 2 points, 51%-75% for 3 points, and 76%-100% for 4 points, while the staining intensity was classified into low, medium, and strong, corresponding to 1, 2, and 3 points, respectively.

The expression levels of S100A10 in both LIHC and normal tissues were assessed using the Wilcoxon rank-sum test. Prognostic analysis was performed using the log-rank test to evaluate survival differences. Additionally, Spearman correlation analysis was conducted to explore the relationship between gene expression and stemness scores. A p-value of less than 0.05 was established as the threshold for statistical significance.

Our investigation commenced with three LIHC samples sourced from the single-cell dataset (GSE112271): GSM3064824, GSM3064820, and GSM3064823. Stringent quality control measures were implemented, requiring each cell to contain a minimum of 200 RNA molecules, a maximum of 2500, and less than 10% mitochondrial RNA ( Figure 1A ). Following this, we employed HARMONY technology to pinpoint highly variable genes from the filtered dataset and conducted bulk deletion analysis based on these feature sets ( Figures 1B–D ). The ANOVA test highlighted ten genes with significant differential expression across the cell samples: TIMP1, IGLL5, CCL21, CXCL10, SPINK1, MT1G, HAMP, SPP1, IGJ, and SAA1 ( Figures 1E, F ). The single-cell analysis categorized the samples into 12 distinct cell populations, including natural killer cells, liver bud hepatocytes, smooth muscle cells, plasma cells, MKI67+ precursor cells, endothelial cells, monocytes, intestinal epithelial cells, cancer stem cells, adventitial cells, dendritic cells, and cancer-associated fibroblasts ( Figures 1G, H ). Notably, functional analysis indicated that stem cell populations are linked to processes such as tumor proliferation, angiogenesis, and epithelial-mesenchymal transition (EMT) ( Figure 1I ).

Initially, we examined expression heatmaps of 16 genes in both TCGA-LIHC samples and normal prostate tissue ( Figure 2A ). We further assessed the relationship between these genes and the clinicopathological features of LIHC patients, visualized through a heat map ( Figure 2B ). The Friends analysis aimed to construct a gene interaction network, utilizing network topology to evaluate gene importance, with MARCKSL1 emerging as a central figure ( Figure 2C ). Univariate COX regression analysis in the TCGA-LIHC dataset demonstrated the prognostic relevance of these 16 genes ( Figure 2D ). We found significant expression differences among SOX4, SH3GL1, RAB11A, FKBP1A, ARL4A, UNC5B, MARCKS, MARCKSL1, LIMS1, STMN1, and LOX across various T phases. Moreover, expression levels of SH3GL1, RAB11A, FKBP1A, and SHC1 differed significantly at varying N stages, and MPP3, SOX4, SH3GL1, RAB11A, FKBP1A, ARL4A, UNC5B, MARCKS, MARCKSL1, LIMS1, STMN1, and LOX showed significant variation across different clinical stages ( Figures 2E–J ). Utilizing the Gene Set Cancer Analysis (GSCA) database, we explored the roles of these 16 genes in LIHC, revealing their involvement in EMT, activation of the cell cycle, and inhibition of estrogen receptor and receptor tyrosine kinase (RTK) pathways ( Figure 2K ). Gene Ontology (GO) analysis indicated that these genes were predominantly associated with calcium-dependent protein binding, membrane microdomain, and pro-B cell differentiation ( Figure 2L ).

To cluster TCGA-LIHC samples, we applied the NMF clustering technique. Co-expression curves were analyzed to determine the optimal classification method for TCGA-LIHC subgroups. The optimal grouping was indicated by the point where the co-expression index exhibited a pronounced decrease. Our results suggested that dividing the samples into two groups was most suitable. Heatmaps illustrating sample divisions into two, three, and four groups are presented ( Figures 3A, B ). We further examined prognostic differences among the groups categorized into two, three, and four clusters. Regardless of the division, patients in cluster 1 consistently demonstrated the poorest prognosis, with p-value analysis confirming that a two-group division was optimal ( Figures 3C–E ). Differential expression of the 16 CSC-related genes across the clusters was also illustrated ( Figures 3F–H ).

Immune checkpoint blockade (ICB) therapy has transformed cancer treatment. In this study, we utilized the TIDE algorithm, which focuses on tumor immune dysfunction and exclusion, to predict the effectiveness of immune checkpoint inhibitors for each TCGA-LIHC sample ( Figure 4A ). TIDE assesses two mechanisms of tumor immune evasion: damage to tumor-infiltrating cytotoxic T lymphocytes (CTLs) and CTL resistance to immunosuppressive factors. Elevated TIDE scores correlate with diminished ICB effectiveness and reduced survival rates following ICB treatment. Upon clustering TCGA-LIHC samples into two groups, we observed differences in ICB response, with patients in cluster C1 showing poorer responses. In three or four cluster categorizations, patients in cluster C2 exhibited better responses to ICB therapy ( Figures 4B, C ). We also employed xCell to analyze immune cell infiltration levels across the TCGA-LIHC samples in different clusters, revealing significant differences in various immune cell types, including CD4+ memory T cells, naive CD8+ T cells, common lymphoid progenitors, M2 macrophages, and plasma B cells ( Figures 4D–G ).

Machine learning is increasingly pivotal in biomedicine, particularly in tumor diagnosis and treatment. By analyzing genomic data from patients, machine learning can facilitate personalized treatment approaches. In this study, we developed a diagnostic model centered on liver hepatocellular carcinoma (LIHC) to aid in the early identification of affected individuals. The model was trained using the TCGA-LIHC dataset and validated with two additional datasets: GSE112790 and GSE102451. Among the 113 tested algorithm combinations, the glmBoost+GBM pair demonstrated the highest efficacy for model construction. For reader convenience, we expanded the prediction results of the top 15 algorithm combinations ( Figure 5A ). The area under the curve (AUC) for the TCGA-LIHC training data was 0.999, while the corresponding AUC values for the validation datasets GSE112790 and GSE102451 were 0.983 and 0.832, respectively. The glmBoost+GBM algorithm identified six key genes: STMN1, SHC1, S100A10, FABP5, ANXA2, and SH3GL1. We subsequently presented the individual predictive values of these genes in TCGA-LIHC, GSE112790, and GSE102451 ( Figures 5B–D ).

Utilizing the XGBoost algorithm, we analyzed the prognostic significance of stem cell-related genes in relation to clinical data from the TCGA-LIHC and IGCG-LIHC datasets to assess their impact on overall survival. Consistently, S100A10, FKBP1A, RAB11A, SOX4, and ARL4A ranked among the top 10 genes in both datasets ( Figures 6A, B ). GOsemSim employs Gene Ontology annotation data to evaluate similarities between gene sets based on shared functional terms. Among these five genes, S100A10 emerged as the most significant ( Figure 6C ). To quantify tumor stemness in LIHC samples, we utilized a logistic regression-based machine learning algorithm (OCLR) described in a Cell article, which computes the stemness index for various samples. Our findings indicated that S100A10, FKBP1A, FABP5, and STMN1 positively correlated with the stemness score of LIHC samples, underscoring the pivotal role of S100A10 in LIHC stemness characteristics ( Figures 6D, E ). Additionally, we analyzed the correlation between S100A10 expression and immune cell infiltration in the LIHC microenvironment, revealing significant differences in microenvironment scores and levels of M2 macrophages, hematopoietic stem cells, endothelial cells, CD4+ Th2 T cells, and NK T cells between high and low expression groups of S100A10 ( Figure 6F ).

The samples within the TCGA-LIHC dataset were divided according to the median expression of S100A10. Those exhibiting expression levels above this median were designated as the high S100A10 expression group, whereas samples with expression below the median were categorized in the low S100A10 expression group ( Figures 7A, B ). To explore the functional implications of S100A10, GO analysis was conducted on differentially expressed genes. S100A10 showed the strongest associations with tubulin binding and cadherin binding in the molecular function (MF) module among the upregulated genes. Additionally, it was linked to the processes of establishing protein localization to the endoplasmic reticulum and nuclear division within the biological process (BP) module. In terms of the cellular component (CC) module, it was connected to cell−substrate junctions and chromosomal regions. On the other hand, among the downregulated genes, S100A10 showed associations with anion transmembrane transporter activity and active transmembrane transporter activity in the MF category, along with involvement in carboxylic acid biosynthesis and organic acid biosynthetic processes in the BP category, and connections to the collagen−containing extracellular matrix and mitochondrial matrix in the CC category, reflecting a very strong relationship. Utilizing KEGG enrichment analysis is an effective approach for unraveling gene functions and providing valuable genomic insights. Importantly, S100A10 demonstrated significant relationships with ribosome and cell cycle pathways among genes that were highly expressed, whereas it was primarily linked to complement and coagulation cascades, as well as the PPAR signaling pathway among genes that were less expressed ( Figures 7C–F ). Subsequently, the ssGSEA algorithm was utilized to systematically compute the enrichment scores for each TCGA-LIHC sample across different pathways, thus establishing a connection between the samples and their corresponding pathways. By evaluating the relationship between gene expression and pathway scores, our goal was to clarify how each gene relates to its associated pathway. The analysis disclosed that S100A10 expression showed a positive correlation with tumor proliferation, G2M checkpoint activity, epithelial-mesenchymal transition (EMT), DNA replication, and DNA repair processes, whereas a negative correlation was found with galactose metabolism ( Figure 7G ). In addition, we investigated the interaction of S100A10 with therapeutic agents for LIHC, and the findings demonstrated that S100A10 has a substantial binding affinity to these drugs, underscoring its potential as a therapeutic target for liver cancer ( Figure 7H ).

In the TCGA-LIHC dataset, we analyzed the correlation between stem cell markers (SOX2, CD44, CD133, and POU5F1) and S100A10. Our findings revealed a significant positive correlation between the expression of S100A10 and CD44, CD133, and POU5F1, while no significant correlation was observed with SOX2 ( Figures 8A–D ). Notably, the expression of S100A10 exhibited the strongest correlation with POU5F1, prompting us to further investigate the relationship between these two proteins. Utilizing The Human Protein Atlas database, we discovered a correlation between the expression levels of S100A10 and POU5F1 at the protein level; however, the limited number of cases precluded statistical significance ( Figure 8E ). To enhance our analysis, we collected 185 liver cancer samples (92 from liver cancer patients and 93 from normal liver tissues) from Nantong Cancer Hospital. This analysis confirmed a positive correlation trend between the expression of S100A10 and POU5F1 at the protein level ( Figures 8F, G ). Lastly, given that POU5F1 is a transcription factor, we examined the potential transcriptional regulatory relationship between S100A10 and POU5F1. Our results indicated that POU5F1 was significantly enriched in the promoter region of S100A10 ( Figure 8H ).

We also analyzed the expression levels of S100A10 in liver cancer samples, revealing that its expression was significantly higher in liver cancer compared to normal liver tissue ( Figures 9A, B ). Additionally, we assessed the predictive value of S100A10 for diagnosing liver cancer using the ROC curve. Our findings indicated that S100A10 possesses a notable predictive value for liver cancer diagnosis (AUC=0.686) ( Figure 9C ). Furthermore, the KM curve analysis demonstrated that patients with elevated S100A10 expression have a poorer prognosis, and S100A10 expression can also serve as a predictor for the 1-year, 2-year, and 3-year prognosis of liver cancer patients ( Figures 9D, E ).