Data for the LIHC samples were obtained from TCGA and the Genotype-Tissue Expression (GTEx, http://www.gtexportal.org/) (20) projects and the fibrosis gene set was obtained from the MSigDB (21). The gene expression matrix was filtered to exclude genes with >50% zero expression values across all samples to ensure data quality and relevance. Additionally, samples were selected from patients with survival times >0 days to focus on clinically meaningful outcomes. The present study employed ConsensusClusterPlus in R (version 4.2.2; R Core Team; http://www.R-project.org/) to classify LIHC patients based on hepatic fibrosis-related gene expression and identified two distinct clusters: Cluster A and Cluster B. The optimal number of clusters was determined using the consensus cumulative distribution function and delta area plot. An unsupervised hierarchical clustering analysis was then conducted to identify genes positively and negatively associated with the clusters, which were visualized through comprehensive heatmaps that incorporated various clinical features. Dimensionality reduction techniques, including uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding, were used to visualize the distribution of gene expression between the two clusters. A differential expression analysis between Cluster A and Cluster B was performed using DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html), where differentially expressed genes were defined as those with an adjusted P-value <0.05 and an absolute log2-fold change >1. The results are illustrated in a volcano plot.

The present study performed Gene Ontology (GO) enrichment analysis using the ClusterProfiler package in R (https://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) to explore the functional differences between Cluster A and Cluster B. This analysis identified enriched biological processes and molecular functions for each cluster. Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) pathway analysis was performed to identify significantly enriched pathways in both clusters. The single-sample gene set enrichment analysis (ssGSEA) algorithm was employed to quantify hallmark pathways in each patient and a differential analysis was conducted to compare these pathways between Cluster A and Cluster B. In addition, ssGSEA was applied to quantify fibrosis and ferroptosis levels in patients and a correlation analysis was performed to explore the relationships between these scores. A differential analysis of the high- and low-fibrosis and ferroptosis groups was performed and Venn diagrams were constructed to identify intersecting genes between these groups.

The present study performed a univariate Cox regression analysis for all intersecting genes to evaluate their prognostic significance. A forest plot was generated to display the top five significant genes, with CALCR showing the smallest P-value. The prognostic role of CALCR was further evaluated across various types of cancer in TCGA dataset, with a focus on overall survival (OS), the disease-free interval (DFI), disease-specific survival (DSS) and the progression-free interval (PFI). Additionally, another forest plot was created to depict the prognostic significance of CALCR based on OS across different types of cancer. In the context of LIHC, a Kaplan-Meier (KM) survival analysis was conducted to evaluate the effects of CALCR on OS, DSS and the PFI. A prognostic calibration curve was plotted to assess the accuracy of the CALCR predictions against the actual outcomes. Finally, a nomogram incorporating CALCR expression and clinical features (T, N and M stages) was constructed to predict patient outcomes.

The present study used the TISCH2 single-cell RNA sequencing database to investigate the single-cell expression profile of CALCR in LIHC. A Pearson's correlation analysis was performed to identify genes that were highly associated with CALCR across different LIHC datasets. Additionally, the correlation of CALCR expression with various cell lines was analyzed within these datasets. The distribution of cell types in each LIHC dataset and the expression pattern of CALCR across these cell types were also examined.

A comprehensive analysis of CALCR expression was conducted across various types of cancer using TCGA data to identify differential expression patterns. The ssGSEA algorithm was employed to quantify hallmark pathway scores for patients across these cancer types, followed by a correlation analysis using Pearson's method to explore the correlation between CALCR expression and these pathways. For LIHC patients, GO and KEGG pathway analyses were conducted to identify significant pathways and biological processes that were differentially enriched between the high- and low-CALCR expression groups. This analysis allowed the further investigation of the functional implications of CALCR expression in various types of cancer and its specific role in LIHC.

The present study performed Pearson's correlation analyses between CALCR expression and RNA modification-related genes in LIHC, focusing on three types of RNA modifications: m1A (10 genes), m5C (13 genes) and m6A (21 genes). Additionally, the associations between CALCR expression and genomic heterogeneity features, including the tumor mutational burden (TMB), microsatellite instability (MSI), neoantigen load (NEO) and homologous recombination deficiency (HRD) were investigated across multiple types of cancer, with a particular focus on LIHC. In LIHC, the correlation between CALCR expression and these genomic features was specifically analyzed. A mutation analysis was conducted to compare CALCR expression in the copy number variation (CNV) loss and CNV gain groups across different cancers, including LIHC. the mutation landscape of CALCR in a pan-cancer context was also examined and differentially mutated genes between the high- and low-CALCR expression groups in LIHC assessed.

The ESTIMATE algorithm (https://bioinformatics.mdanderson.org/estimate/) was applied to calculate stromal, immune and ESTIMATE scores for LIHC patients. Correlation analyses were performed using Pearson's method to investigate the relationships between these scores and CALCR expression, providing insights into the interaction of the tumor microenvironment with CALCR. Additionally, a differential expression analysis of immune checkpoint genes between the high- and low-CALCR expression groups was conducted, with a focus on both inhibitory and stimulatory checkpoint genes. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used and a SubMap analysis performed to compare the potential responses of the high and low CALCR expression groups to immune checkpoint blockade therapy (22,23). These findings provided a comprehensive view of the role of CALCR in modulating the immune response and its potential as a predictor of immunotherapy efficacy.

The present study conducted a differential expression analysis of CALCR using data from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) and Human Protein Atlas (HPA, http://www.proteinatlas.org/) databases to validate the bioinformatics findings. Data for CALCR expression were retrieved from the GEO database (GSE142987; public on September 30, 2020) to assess its differential expression in LIHC and normal tissues. The analysis was conducted using a predefined dataset and statistical significance was determined using a t-test and P<0.05 was considered to indicate a statistically significant difference. The expression levels of CALCR were further analyzed across different tumor stages. Immunohistochemical data from the HPA were analyzed to assess CALCR expression in both normal and LIHC tissues. Tissue sections from a normal liver (female, age 73) and a liver with LIHC (male, age 73) were examined for CALCR staining intensity. The presence and localization of CALCR expression were evaluated in both the cytoplasmic and membranous regions, with the staining intensity categorized as negative, low, moderate, or high. These analyses provide a comprehensive evaluation of CALCR expression across normal and cancerous tissues, reinforcing its potential role in tumorigenesis.

Dulbecco's modified Eagle's medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.) was used to culture the liver cancer cells in a moist environment at 37°C and 5% CO₂, including HepG-2 and HuH-7 cell lines, which were obtained from the American Type Culture Collection. Lipofectamine^® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) was used to transfect synthesized small interfering (si)RNA (Guangzhou RiboBio Co., Ltd.). The siRNA was dissolved in a solution with a concentration of 20 µM, and the transfection efficiency was detected after transfection for 24 h at 37°C. The cells were successfully transfected and continued to be cultured for 24 h before other experiments were performed. The effectiveness of silencing was evaluated via reverse transcription-quantitative (RT-q) PCR, with the reference gene indicated as GAPDH. The CALCR siRNA sequences are shown in Table SI.

Total RNA was extracted from liver cancer cells (1×10⁷ cells) using TRIzol^® reagent (Thermo Fisher Scientific, Inc.) according to the manufacturer's protocols. The purity and concentration of total RNA were determined by measuring the absorbance values at 260 and 280 nm. Then, cDNA was synthesized using high-capacity cDNA reverse transcription kits (Vazyme Biotech Co., Ltd.) following the manufacturer's instructions. Finally, SYBR Green PCR Master Mix (Vazyme Biotech Co., Ltd.) was used to perform RT-qPCR according to the manufacturer's protocols, with the reference gene indicated as GAPDH. The CALCR primers used are listed in Table SII. The thermocycling conditions were: Initial denaturation at 95°C for 10 min, followed by 40 cycles of 95°C for 15 seconds and 60°C for 1 min. The method of quantification was based on the comparative Cq method (24). The experiments were conducted with four biological replicates and three technical replicates.

After culture for 24, 48, or 72 h, 10 µl of the Cell Counting Kit-8 (CCK-8) reagent was added to cells in the 96-well plates and incubated for an additional 4 h in the dark. The absorbance values were measured using an automated microplate reader (Synergy4; BioTek, USA) at 450 nm.

For the migration assay, 5×10⁵ cells suspended in serum-free medium were seeded in the upper chamber and medium containing 20% FBS was added to the lower chamber. The membranes in the chambers were fixed with 4% paraformaldehyde at 37°C for 30 min and stained with 0.1% crystal violet at 37°C after a 24 h incubation.

For the wound healing assay, cells were seeded into a 6-well plate and scratched with a 1 ml pipette tip after the confluence reached 90%. The diameter of the wounds was assessed after 24 h of treatment with 1 µg/ml mitomycin C at 37°C in serum-free media using an inverted microscope (ICX41; SDPTOP). ImageJ software (Version 1.54f 29 June 2023; National Institutes of Health) was used to measure the wound widths and wound closure rates, with images of the wounds captured at 0 and 24 h after scratching.

An Annexin V-fluorescein isothiocyanate (FITC)/propidium iodide (PI) apoptosis kit (Multi Sciences) was used to detect cell apoptosis. The cells were collected and rinsed twice with cold phosphate-buffered saline, followed by staining with 5 µl of Annexin V-FITC for 15 min and 5 µl of PI for 5 min at 37°C. Cell apoptosis was identified using a CytoFLEX-3 flow cytometer (Beckman Coulter, USA) and the collected data were analyzed via FlowJo software (v10.8.1; FlowJo LLC). The apoptotic rate was calculated as the percentage of early + late apoptotic cells.

All statistical analyses were conducted using R software (version 4.2.2). The data from the TCGA database were analyzed using the Mann-Whitney U test, while the data from the GEO database were analyzed using Student's t-test. For comparisons involving more than two groups, one-way ANOVA was applied, followed by the Bonferroni's post hoc test to correct for multiple testing. All P-values were adjusted using the FDR correction to account for multiple comparisons. P<0.05 was considered to indicate a statistically significant difference. The Kruskal-Wallis test followed by Dunn's test was employed to evaluate the association between clinical parameters and CALCR expression. Survival analyses, including OS, DSS, DFI and PFI, were conducted using Kaplan-Meier curves and significance was determined using the log-rank test. The prognostic value of CALCR in patients with LIHC was evaluated via univariate Cox regression analyses. The patients were split into high and low expression groups by the median. Finally, the correlation between CALCR expression and clinical characteristics was evaluated using Pearson's correlation coefficient. In cell experiments, the two-group analysis was conducted via the one-way analysis of variance with the Student's t-test. All experiments were independently repeated at least three times and data are presented as the mean ± standard deviation.

The consensus clustering analysis effectively stratified the LIHC patients into two distinct clusters: Cluster 1 and Cluster 2 (Fig. 1A). Unsupervised clustering and Pearson's correlation analyses revealed a clear segregation of genes, with the results supported by the integration of various clinical features (Fig. 1B). Visualization in UMAP and tSNE plots revealed a pronounced separation between the clusters, indicating significant differences in gene expression profiles (Fig. 1C and D). The Kaplan-Meier survival analysis indicated that patients in Cluster A had an improved overall survival prognosis compared with those in Cluster B (P=0.0034; Fig. 1E). A differential expression analysis using DESeq2 identified differentially expressed genes (DEGs) between Cluster A and Cluster B. These DEGs were defined as those with an adjusted P-value <0.05 and an absolute log2 fold change >1. The volcano plot illustrates the molecular heterogeneity between the two clusters (Fig. 1F). These findings underscore the distinct molecular landscapes of the identified clusters, providing a foundation for further exploration of the biological and clinical implications of hepatic fibrosis in LIHC patients.

GO enrichment analysis revealed that Cluster A was significantly enriched in pathways related to iron ion binding and transition metal ion binding (Fig. 2A), whereas Cluster B was enriched in pathways associated with ion channel activity and metal ion transmembrane transporter activity (Fig. 2B). KEGG pathway analysis identified distinct enriched pathways in both clusters (Fig. 2C and D). The hallmark pathway analysis using ssGSEA quantified pathway scores for each patient and the differential analysis highlighted significant differences between Cluster A and Cluster B (Fig. 2E). The correlation analysis of the ssGSEA scores revealed a positive correlation between fibrosis levels and ferroptosis levels in patients (P<0.001; Fig. 2F). Furthermore, the differential expression analysis of the high- and low-fibrosis groups and the high- and low-ferroptosis groups revealed intersecting genes, which were visualized in Venn diagrams for both the high-fibrosis/high-ferroptosis (Fig. 2G) and low-fibrosis/low-ferroptosis groups (Fig. 2H). These comprehensive analyses underscored the intricate functional differences and correlations within fibrosis-related clusters in LIHC patients, providing valuable insights into the molecular foundations of these conditions.

The univariate Cox regression analysis identified five significant genes from the intersecting gene set, with CALCR showing the most significant prognostic value (P=0.002; Fig. 3A). The prognostic role of CALCR was evaluated in TCGA pan-cancer dataset, indicating its effects on OS, DFI, DSS and PFI across multiple types of cancer (Fig. 3B). The forest plot for OS across different cancers highlighted the significant prognostic role of CALCR (P<0.001; Fig. 3C). The KM survival analysis showed that high CALCR expression was associated with shortened OS (P=0.002) and DSS (P=0.013) in LIHC patients but not with the PFI (P=0.457; Fig. 3D). The prognostic calibration curve for CALCR showed good alignment with the ideal curve, indicating reliable predictive performance (Fig. 3E). A nomogram incorporating CALCR expression and clinical features (T, N and M stages) was constructed to predict patient outcomes, providing a practical tool for clinical prognostication (Fig. 3F). These findings underscored the significant prognostic value of CALCR in both pan-cancer and LIHC contexts, highlighting its potential as a biomarker for patient stratification and outcome predictions.

The single-cell analysis of the TISCH2 database revealed that multiple genes were strongly associated with CALCR across different LIHC datasets (Fig. 4A). While the correlation analysis using Pearson's method did not reveal significant associations between CALCR expression and various cell lines in these datasets (Fig. 4B), the distribution of cell types in the LIHC datasets and the expression pattern of CALCR within these cell types were visualized. This visualization revealed distinct expression profiles across different cell populations (Fig. 4C). These findings provided insights into the single-cell expression landscape of CALCR in LIHC and its potential interactions with other genes and cell types.

The pan-cancer differential expression analysis revealed that CALCR was significantly differentially expressed across a variety of cancer types, including LIHC (P<0.001; Fig. 5A). In LIHC, a correlation analysis based on ssGSEA pathway scores revealed that CALCR expression was positively associated with pathways such as the G2/M checkpoint and E2F targets, which are involved in cell cycle regulation. By contrast, CALCR expression was negatively associated with pathways involved in bile acid metabolism, fatty acid metabolism and xenobiotic metabolism, which are essential for metabolic processes (Fig. 5B). Furthermore, GO and KEGG pathway analyses of LIHC patients revealed significant differences in pathways and biological processes between the high- and low-CALCR expression groups, emphasizing the multifaceted roles of CALCR in different cellular contexts (Fig. 5C and D). These results provided a comprehensive overview of the involvement of CALCR in various types of cancer, with specific functional implications in LIHC.

The Pearson's correlation analysis in LIHC revealed significant associations between CALCR expression and 44 marker genes involved in RNA modifications, including m1A, m5C and m6A modifications (Fig. 6A). The genomic heterogeneity analysis demonstrated that CALCR expression was positively associated with the TMB, MSI, NEO and HRD across various types of cancer, with particularly strong correlations observed in LIHC (Fig. 6B). The mutation analysis revealed similar CALCR expression levels between the CNV loss and CNV gain groups of LIHC patients, although differences were noted in patients with other cancers (Fig. 6C). The mutation landscape analysis of CALCR across various types of cancer identified specific mutations associated with CALCR expression (Fig. 6D). In LIHC, differentially mutated genes between the high- and low-CALCR expression groups were visualized in a waterfall plot, highlighting distinct mutation profiles (Fig. 6E). These analyses reveal complex interactions between CALCR expression and various genomic features, providing insights into its role in cancer biology.

The correlation analysis using the ESTIMATE algorithm revealed a negative correlation between ESTIMATE scores and CALCR expression, suggesting a more complex immune cell infiltration environment in patients with low CALCR expression Specifically, the correlation analysis showed significant associations, with StromalScore (P=0.018), ImmuneScore (P<0.001) and ESTIMATEScore (P<0.001) all demonstrating negative correlations with CALCR expression (Fig. 7A-C). The differential expression analysis of immune checkpoint genes showed significant differences between the high and low CALCR expression groups, with notable variations in both inhibitory (Fig. 7D) and stimulatory (Fig. 7E) checkpoint genes. The TIDE analysis revealed no significant difference in overall TIDE scores between the high- and low-CALCR expression groups; however, exclusion levels were significantly lower in the low-CALCR expression group (Fig. 7F). The SubMap analysis indicated that patients with low CALCR expression might exhibit an improved response to immunotherapy (Fig. 7G). These results elucidate the complex interplay between CALCR expression and the immune microenvironment in LIHC, highlighting potential implications for immunotherapy efficacy.

The present study conducted an experimental validation of CALCR expression using data from the GEO and HPA databases to further corroborate the bioinformatics findings. These analyses were designed to assess the differential expression of CALCR in cancerous and normal tissues, thereby reinforcing the role of this gene in tumorigenesis.

First, an analysis of the GSE142987 dataset, which included 35 liver cancer patients and 30 healthy controls, revealed a statistically significant difference in CALCR expression between tumor and normal tissues. Specifically, CALCR expression was markedly elevated in tumor samples compared with normal controls (P=0.033). These findings suggested the potential upregulation of CALCR expression under cancerous conditions, reinforcing its possible involvement in cancer development and progression. In the analysis of CALCR expression across different tumor stages using the GSE142987 dataset, significant differences were observed between specific stages. Stage: n/a showed a statistically significant difference in CALCR expression compared with Stage: A (P<0.01), as did Stage: A when compared with Stage: B/C (P<0.001). However, CALCR expression in Stage: 0 did not exhibit a significant difference compared with other stages, with all comparisons yielding P-values >0.05. These results highlighted the differential expression of CALCR in later tumor stages while suggesting minimal variation in early-stage tumors (Fig. S2A and B).

Further validation was conducted using the HPA database, which focuses on the immunohistochemical staining of CALCR in both normal and LIHC tissues. In normal liver tissue (female, age 73), CALCR expression was not detected and staining with the antibody was not observed. By contrast, in LIHC tissue (male, age 73), CALCR expression was detected at low levels, with a moderate staining intensity observed in <25% of the tissue that was localized primarily in the cytoplasmic and membranous regions (Fig. S2C-D). These results, despite the low expression in cancerous tissues, further supported the differential expression of CALCR between normal and cancerous tissues, highlighting its potential relevance in the oncogenic context.

The combined findings from the GEO and HPA analyses provide strong evidence for the differential expression of CALCR between normal and cancerous tissues, underscoring its importance in liver cancer and possibly other cancer types.

In HepG-2 and HuH-7 cells, si-CALCR-2 effectively reduced CALCR expression (Fig. 8A and B). Compared with the si-NC group, the CCK-8 analysis revealed that CALCR knockdown significantly affected the proliferation of both HepG-2 and HuH-7 cells (Fig. 8C and D). Additionally, CALCR depletion substantially inhibited the migration of HepG-2 and HuH-7 cells, as demonstrated by Transwell migration (magnification, ×200; scale bar, 100 µm) and wound healing (magnification, ×100; scale bar, 120 µm) assays (Fig. 8E-H). According to the flow cytometry findings, CALCR markedly reduced apoptosis in LIHC cells (Fig. 8I).