This study utilized data from TCGA (https://portal.gdc.cancer.gov) (50 samples of healthy liver and 374 samples of HCC) and GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14520) (242 HCC samples) databases to collect clinical information and transcription profiling data of patients with HCC. Gene IDs were converted into corresponding gene symbols using the human gene annotation file, with mean values used for multiple probes targeting the same gene ID. The TCGA dataset was used as the training set, while the GEO dataset served as the testing set. A total of 224 BMRGs previously identified in other studies were retrieved ( Supplementary Table 1 ) (18), and their expression matrix was extracted for subsequent analysis after data matching, filtering, and correction.

The “Limma” R package was utilized for differential expression analysis of BMRGs in normal and cancer tissue samples, with statistical significance defined as FDR< 0.05 and LogFC< 0.585, resulting in the identification of 585 genes exhibiting significant differential expression. Subsequently, the differentially expressed genes (DEGs) underwent GO and KEGG pathway enrichment analysis using the “Clusterprofiler” R package to elucidate their biological characteristics and cellular functional pathways. Statistical significance was determined based on P-values and adjusted P-values< 0.05. Finally, visualization of the enrichment analysis results was performed using the “ggplot2” and “goplot” R packages. Gene mutation data were obtained from the TCGA database to calculate tumor mutation burden (TMB) in HCC patients, while exploration of DEG mutation frequency was conducted using the “Maftools” R package.

First, the prognostic outcomes were integrated with the expression levels of DEGs in each sample. In the training set, univariate Cox regression analysis was employed to identify genes associated with prognosis among DEGs. Genes with P-values< 0.05 were selected for analyzing the association between gene mutation frequency and mutated genes in HCC samples from the training set using the “Maftools” R package. Subsequently, a Least Absolute Shrinkage Selection Operator (LASSO) Cox regression analysis was performed using the “glmnet” R package to further refine prognostic-related genes and develop a prognostic risk score model for predicting overall survival (OS) in HCC samples. The risk score for each sample was calculated using the following formula.

The “Coef” represents non-zero regression coefficients obtained through LASSO Cox regression analysis ( Supplementary Table 2 ), while “ExpGene” denotes the expression value of a gene derived from a prognostic risk score model. All samples were stratified into high-risk and low-risk groups based on the median risk score, and Kaplan-Meier analysis was conducted using the log-rank test to compare overall survival differences between these two groups. Time-dependent receiver operating characteristic (ROC) curves were generated using the “Survival ROC” package in R to evaluate the predictive accuracy of prognostic risk scoring patterns. Finally, the validity and applicability of the prognostic risk scoring model were further assessed in an independent test set.

To comprehensively investigate the distinct disparities between the high-risk and low-risk score groups, we successfully employed the “Limma” package in R to conduct principal component analysis (PCA) on the gene expression profiles before and after implementing the prognostic risk score model within our training dataset. Initially, PCA was performed on all DEGs associated with the basement membrane, followed by an examination of gene expression profiles derived from the prognostic risk score model using PCA. Ultimately, these PCA results were visually presented on a two-dimensional graph utilizing the “ggplot2” package.

The risk scores for each sample were integrated with the corresponding clinical information, encompassing gender, age, pathological stage, and TNM stage. The relationship between the risk score and clinical characteristics was investigated using the “limma” R package. Based on these clinical characteristics, the samples were stratified into two groups to facilitate comparison of differences in risk scores. Statistical significance was determined by considering P-values< 0.05.

The overall survival of individual patients with HCC was further investigated by developing nomograms using the “rms” package in R, incorporating age, gender, pathological stage, and prognostic risk score models. The predictive ability of these nomograms was evaluated through ROC curves and calibration plots.

Gene Set Enrichment Analysis (GSEA) was conducted on the “c2.cp.kegg.v7.4.symbols” gene sets obtained from the molecular signature database using the “org.Hs.eg.db” and “clusterProfiler” packages to assess differences in biological processes between high-risk and low-risk groups (https://www.gsea-msigdb.org/gsea/msigdb/). Additionally, single-sample gene set enrichment analysis (ssGSEA) was performed using the “GSVA” and “GSEABase” R packages to visualize infiltration fractions of 16 immune cells and activity levels of 13 immune-related pathways as ssGSEA scores. To predict the potential impact of immune checkpoint blockade, we also examined expression levels of PD-1, PD-L1, CTLA4, PD-L2, Havcr2, and TIGIT (19). Furthermore, we utilized the “PRRophetic” R package to predict semi-inhibitory concentrations of Sorafenib, Sunitinib, Erlotinib, and Gemcitabine in each sample as indicators for their effectiveness in inhibiting specific biological or biochemical functions (20). Lastly, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) online database (http://TIDE.dfci.harvard.edu/) to forecast immunotherapy efficacy in both high-risk and low-risk populations; statistical significance was considered when P-value< 0.05 (21).

Firstly, the transcription profiling data of high-risk and low-risk scoring groups were compared using the ‘limma’ R package. Differentially expressed genes were identified based on a false discovery rate<0.05. The differentially expressed genes were then analyzed using the STRING online database (version: 11 5; https://cn.string-db.org/) to generate a PPI network with an interaction score >0.70 (high confidence). Subsequently, Cytoscape software (version: 3.9.1) was utilized for further processing and visualization of the PPI network data. The top 10 most central genes were screened using the cytoHubba plugin (version: 0.1). Next, all samples were stratified into high-expression and low-expression groups using the ‘cut off’ R package. Kaplan-Meier analysis was performed to assess survival differences between these two groups. Finally, we investigated immune cell infiltration patterns and clinical features associated with prognosis-related central genes.

The Wilcoxon rank-sum test was employed to compare differences between the two groups. Kaplan-Meier analysis was utilized to evaluate survival disparities between high-risk and low-risk score groups. Multivariable Cox regression analysis was conducted to identify independent predictors of OS in HCC. ROC curves were plotted to assess the predictive validity of prognostic risk scoring models and nomograms. All statistical analyses were performed using R 4.0.2, with a significance level set at P< 0.05.

Additional materials and methods details are provided in the online Supplementary Material .

We present the study’s flow chart in Supplementary Figure 1 . A comprehensive analysis of the expression levels of 224 BMRGs in tumor and normal samples from the TCGA database, we identified 121 DEGs in the TCGA HCC cohort (P< 0.05, FDR< 0.585). Among these, 113 genes were upregulated, and 8 genes were downregulated in tumor samples. Figures 1A and B illustrate the heat map and volcano plot of DEGs between tumor and normal samples, respectively.

The biological functions of the differentially expressed BMRGs were further determined through GO annotation and KEGG pathway enrichment analyses. The identified BMRGs were primarily associated with extracellular matrix organization, extracellular structure organization, and cell-substrate adhesion in terms of biological processes (BP). In relation to cellular components (CC), the BMRGs were mainly enriched in collagen-containing extracellular matrix, basement membrane, and endoplasmic reticulum lumen. Regarding molecular functions (MF), the BMRGs predominantly exhibited roles as extracellular matrix structure constituents and cell adhesion molecule binders ( Figure 1C ). The top 80 significantly enriched genes and pathways are displayed in Figure 1D . KEGG pathway enrichment analysis revealed significant enrichments for ECM-receptor interaction, focal adhesion, human papillomavirus infection, and PI3K-Akt signaling pathway ( Figure 1E ). Figure 1F presents the top 45 significantly enriched genes and pathways. Collectively, these findings indicate a substantial correlation between BMRGs and extracellular matrix-associated functions and pathways in HCC progression.

The TCGA-HCC cohort was utilized as the training set in this study, leading to the identification of 121 BMRGs. Univariate Cox analysis revealed 31 BMRGs significantly associated with OS ( Figure 2A ). Analysis of the somatic mutation spectrum in these prognostic genes demonstrated a mutation frequency of 24.53% in the HCC samples (n=371), with a total of 91 cases exhibiting mutations ( Figure 2B ). ROBO1 exhibited the highest mutation frequency, followed by ADAMTS9, ITGB5, SLIT1, Lamb1, ADAMTS7, CSPG4, POSTN, ITGAM, ADAM17, ITGA3, LAMC1, MEP1A, ITGA2, ITGAV, NID1, LAMA4, MMP14, PTPRF, ROBO3, SMC3 and ITGB1, while no mutations were observed in other genes. Further analysis revealed positive associations between Robo1 and ITGAM, ADAMTS9 and CSPG4, ITGA2 and LOXL2, NID1 and CD151 ( Figure 2C ). These 31 genes were included in a LASSO Cox regression analysis based on the TCGA-HCC cohort. The analysis identified six genes, namely CD151, CTSA, MMP1, ROBO3, ADAMTS5, and MEP1A, which were used to construct the risk score model for BMRGs ( Figures 2D, E ). The risk score for each sample was calculated using the following formula: risk score = (0.00102 × expression of CD151) + (0. 192 × expression of CTSA) + (0. 165 × expression of MMP1) + (0.0739 × expression of ROBO3) + (0.302 × expression of ADAMTS5) + (0.0203 × expression of MEP1A).

The established risk score model based on BMRGs was employed to effectively stratify the HCC cohorts into two distinct subgroups: high-risk and low-risk groups ( Figures 3A, B ). In this risk model, TCGA-HCC patients were classified into a high-risk group (n = 185) and a low-risk group (n = 185) using the median risk score as the cutoff, while GEO-HCC patients were divided into a high-risk group (n = 111) and a low-risk group (n = 110). Notably, in both the training and test sets, the low-risk group exhibited significantly improved clinical outcomes compared to the high-risk group (P-value< 0.05) ( Figures 3C, D ). The distribution of risk scores in relation to age, gender, pathological stage, and TNM stage of HCC was analyzed. Although there were no significant differences observed in the association between risk scores and age, gender, or M stage, higher risk scores were associated with elevated T, N, and pathological stages ( Supplementary Figure 2 ). Univariate prognostic Cox analysis revealed that both pathological stage and risk score independently served as prognostic factors. This finding was further confirmed by multivariate Cox analysis (P-values< 0.05) ( Figures 3E, F ). The model’s reliability was assessed using ROC curves, yielding area under the curve (AUC) values of 0.773, 0.695, and 0.643 for years 1, 3, and 5, respectively ( Figure 3G ). Additionally, the risk model demonstrated an AUC of 0.643 ( Figure 3H ). These findings provide evidence supporting the reliability of the basement membrane risk model in predicting HCC progression.

To accurately predict the probability of OS, we developed a nomogram that integrated the BMRGs score with other clinicopathological features, including gender, age, and pathological stage. The predictive performance of the nomogram for 1-year, 3-year, and 5-year OS in patients with HCC was assessed ( Figure 4A ). The calibration curve demonstrated close alignment with the ideal curve, indicating high predictive accuracy ( Figure 4B ). Furthermore, ROC curve analysis was conducted to validate the utility of the nomogram and calculate the AUC for risk assessment as well as for age, gender, and pathological stage prediction ( Figure 4C ). Notably, our nomogram exhibited the highest AUC value among all factors evaluated, indicating its superior prognostic accuracy. Univariate and multivariate Cox analyses performed on data from the training set confirmed that our nomogram served as an independent prognostic factor (P-value< 0.05) ( Figures 4D, E ). Collectively, these findings provide comprehensive validation of both the clinical applicability and predictive power associated with our BMRGs-based prognostic model.

To comprehensively analyze TME, we utilized CIBERSORTx to calculate the infiltration degree of immune cell profiles. Our findings revealed a significant increase in the abundance of activated dendritic cells (aDCs), immature dendritic cells (iDCs), macrophages, T helper 1 (Th1) cells, T helper 2 (Th2) cells, and regulatory T cells (Tregs) in the BMRGs high-risk group. Conversely, the BMRGs low-risk group exhibited a significantly increased abundance of natural killer (NK) cells ( Figure 5A ). Immune function analysis demonstrated that APC co-stimulation, CCR signaling pathway activation, checkpoint regulation, HLA expression modulation, MHC class I antigen presentation enhancement and parainflammation-related functions were highly active in the high-risk group. In contrast, Type II IFN-Response immune-related functions were prominently activated in the low-risk group ( Figure 5B ). Furthermore, our results indicated that the high-risk group had lower Tumor Immune Dysfunction and Exclusion (TIDE) scores compared to the low-risk group ( Figure 5C ), suggesting enhanced effectiveness of immunotherapy within this subgroup. Risk score and drug sensitivity analysis revealed higher drug sensitivity to Sorafenib, Sunitinib, and Gemcitabine among patients classified into the high-risk group as opposed to those assigned to the low-risk category. However, Erlotinib displayed reduced drug sensitivity within the high-risk group ( Figures 5D-G ). Given its significance in checkpoint inhibitor immunotherapy approaches, we also investigated associations between risk scores and key immune checkpoints. Our findings indicated elevated levels of PD-1, PD-L1, CTLA4, PD-L2, HAVCR2, and TIGIT expression within individuals classified into the high-risk group ( Figure 5H ), thereby indicating an existence of immunosuppression within this subgroup. In conclusion, the quantification of BMRGs risk score represents a novel and reliable biomarker for evaluating prognosis associated with immunotherapy.

The interactions between DEGs in high-risk and low-risk groups were analyzed using the STRING online database, revealing a complex PPI network ( Supplementary Figure 3A ). The PPI network data were further processed and visualized using Cytoscape software, highlighting the intricate interactions among DEGs ( Supplementary Figure 3B ). Hub genes were identified using the cytoHubba plugin in Cytoscape, resulting in the selection of ten key genes within the network, namely MUC1, MUC6, MUC5AC, ENO2, PKM, CXCL8, ITGA3, COL3A1, GCNT1, and LIF ( Figure 6A ). These hub genes were sequenced using the degree method, as detailed in Supplementary Table 2 . Differential analysis revealed that PKM and ITGA3 are significantly overexpressed in HCC patients compared to normal samples (P< 0.05) ( Figures 6B, C ). Moreover, both PKM and ITGA3 showed a significant association with survival outcomes, where higher gene expression correlated with worse prognosis (P< 0.05) ( Figures 6D, E ). Furthermore, PKM expression was found to increase with tumor stage progression, indicating its potential involvement in the advancement of HCC ( Figure 6F ). In contrast, ITGA3 expression did not show a significant correlation with tumor stage, although a decreasing trend in advanced stages was observed ( Figure 6G ). These results suggest that PKM and ITGA3 could serve as potential biomarkers and therapeutic targets for HCC.

To further validate the association of PKM2 and ITGA3 with HCC at the clinical level, we first examined the mRNA levels of PKM2 and ITGA3 in 165 paired HCC and corresponding adjacent nontumor specimens by quantitative RT-PCR (RT-qPCR). PKM2 and ITGA3 mRNA levels were significantly elevated in HCC specimens compared to adjacent nontumor specimens ( Figures 7A, B ).

We then assessed the protein levels of PKM2 and ITGA3 in paired HCC specimens. Western blot assay showed that PKM2 and ITGA3 were significantly elevated in HCC specimens compared to adjacent nontumor specimens ( Figures 7C, D ).

To explore the clinical significance of PKM2 and ITGA3, HCC patients were categorized into two groups based on RT-qPCR, negative or positive of PKM2/ITGA3. Notably, Kaplan-Meier analysis revealed that higher PKM2 or ITGA3 levels were associated with worse OS and worse Progression-Free Survival (PFS) in HCC patients ( Figures 7E, F ). Clinicopathological characterization indicated that PKM2 or ITGA3 overexpression was associated with age, serum AFP, TNM stage, BCLC stage, tumor number, tumor size, and microvascular invasion ( Table 1 ). Multivariate analysis demonstrated that PKM2 and ITGA3 are significant and independent predictors of OS and PFS in HCC patients ( Supplementary Tables 3 , 4 ). In conclusion, our data suggest that regulation of PKM2 and ITGA3 may predict poor prognosis in HCC patients.

To investigate the biological functions associated with BMRGs, we employed GSEA to analyze the TCGA cohort. We identified 37 pathways that were significantly enriched in the high-risk group and 8 pathways in the low-risk group (P< 0.05). Notably, the top five enrichment pathways in the high-risk group included cell adhesion molecules (CAMs), cytokine-cytokine receptor interaction, ECM-receptor interaction, focal adhesion, and neuroactive ligand-receptor interaction ( Supplementary Figure 4A ). Conversely, the top five enrichment pathways in the low-risk group comprised beta-alanine metabolism, fatty acid metabolism, linoleic acid metabolism, primary bile acid biosynthesis and retinol metabolism ( Supplementary Figure 4B ). These findings underscore the distinct biological processes that are activated in high-risk and low-risk groups, providing insights into the molecular mechanisms underlying the progression and treatment response of HCC.