Identification of a novel matrix metalloproteinases-related prognostic signature in hepatocellular carcinoma

Datasets

The database Mammalian Degradome (http://degradome.uniovi.es/dindex.html) was used to look for MMPs-related genes. RNA sequencing profiles and patient information for 374 HCC samples were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). In addition, 231 HCC patients from the International Cancer Genome Consortium for the Study of Liver Cancer in Japan (ICGC-LIRI-JP) dataset and their corresponding clinical characteristics were obtained from the ICGC database (https://dcc.icgc.org/). Furthermore, a single-cell sequencing set of HCC (GSE189903) was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The acquisition of all publicly available data adheres to the publication and database access policies for the mentioned databases. A flowchart for bioinformatics analysis of publicly available datasets from the TCGA, ICGC, and GEO databases is shown in Figure 1.

Figure 1. Flowchart for bioinformatics analysis in this study.

Identification of prognostic MMPs-related DEGs

A univariate Cox regression analysis was performed to assess the predictive potential of MMPs-related genes and identify the association between MMPs-related differentially expressed genes (DEGs) and survival status in the TCGA-LIHC cohort.

Consensus clustering of prognostic MMPs-related DEGs

Clustering analysis of the expression patterns of survival-related MMP genes to explore the association between their expression and survival status was performed using the “ConsensuClusterPlus” package. One hundred replicates were performed with a pltem value of 0.8 to verify the stability of the subtypes. The clustering effect was best when the k-value was set to 2. Kaplan-Meier analysis was performed using the “survminer” package to estimate the differences in overall survival (OS) and progression-free survival (PFS) between different clusters. A log-rank test was performed to assess differences in survival. The above survival-related genes between the two clusters were analyzed using the “limma” package in R.

Construction and validation of MMPs-related prognostic signature

For modeling purposes, unnecessary genes were removed by using LASSO Cox regression. The risk score was then calculated using a formula derived from the correlation coefficient and gene expression value obtained from the multifactorial Cox regression analysis. The formula was then applied to determine the scores, which were then used to categorize patients based on the median into high-risk and low-risk groups.

Subsequently, the TCGA cohort was assigned as a training cohort and a testing cohort, with the testing cohort serving as an internal validation cohort and the ICGC as an external validation cohort, both used to validate the predictive model.

To investigate the gene expression distribution in the predictive model, two separate analyses, namely t-distributed stochastic neighbor embedding (t-SNE) and principal component analysis (PCA), were undertaken to assess positioning for two distinct risk groups. We evaluated the OS between training and the testing cohorts by employing Kaplan-Meier analysis and visualized the results using the “survminer” package. The AUC was calculated at 1-, 3-, and 5-year intervals using the “time ROC” package to assess the accuracy of the predictive model.

Risk score independent prognostic analysis

Univariate and multivariate Cox regression analyses assessed the relationship between risk factors and prognosis. Forest maps were plotted to show the independent predictive value of the risk score. A nomogram was created using the risk score and other clinical indicators to predict 1-, 3- and 5-year OS of HCC patients. To further evaluate the discriminative and predictive ability of the nomogram, we also calculated the concordance index (C-index). The C index values varied from 0.5 to 1.0, whereby a higher C index value indicates a higher discriminative ability of the prediction model [14].

Analysis of functional enrichment

The analyses of Gene Ontology (GO) and the Kyoto Encyclopedia of Gene and Genomes (KEGG) were accomplished using the “clusterProfiler” package. Furthermore, we used Gene Set Enrichment Analysis (GSEA) to determine the significance of the MMPs-related prognostic signature in differentiating between the low and high-risk groups.

Evaluation of the infiltration of immune cells

We employed single-sample Gene Set Enrichment Analysis (ssGSEA) as well as Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) to assess the magnitude with which the infiltration of immune cells was observed in different clusters and groups of the immune analysis. Gene set for ssGSEA was acquired via previous study [15]. The number of iterations of CIBERSORT was set to 1000. We excluded samples for which the P-value was less than 0.05. Additionally, the correlation between the immune cell profiles and the MMPs-related prognostic signature was examined using Spearman’s correlation. Based on the TCGA-LIHC profile of solid tumor expression, six immune subtypes (immune C1-C6) were identified, including wound healing, IFN-γ dominant, inflammatory, lymphocyte depleted, immunologically quiet, and TGF-beta dominant [16]. A two-way analysis of variance (ANOVA) was used to examine the correlation between the risk score and immune subtypes.

Analysis of drug sensitivity

Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) is an accessible dataset that provides comprehensive data regarding the sensitivity of cancer cells to drugs as well as drug-response molecular markers [17]. By employing the package of “oncopredict”, we distinguished the responsiveness to different groups of medications, enabling the assessment of 198 therapeutic compounds. The sensitivity scores were then evaluated to determine the predicted IC₅₀ of all drugs in HCC patients.

Profiling of RNA-seq at the single-cell level

The R package “Seurat” was used for unsupervised clustering of single cells from cancer or normal samples, and “Seurat” for annotation. Two clustering methods for dimensionality reduction, tSNE and Uniform Manifold Approximation and Projection (UMAP) were used in this study. The expression of the MMPs-related signature was visualized using the “VlnPlot” function of the “Seurat” R package.

Statistics analysis

All statistical analyses were performed using R software (version 4.2.3). Continuous variables were tested using the student t-test, while categorical variables were tested using the chi-squared test. A P-value < 0.05 was considered significant.

Identification of prognostic MMPs-related DEGs

This study obtained 203 well-defined MMPs-related genes (Supplementary Table 1). A total of 13 prognostic MMPs-related DEGs significantly correlated with OS were identified from TCGA-LIHC cohort by univariate Cox regression analysis based on the criterion of FDR<0.005, including RNPEPL1, MMP1, ADAM9, ADAM15, ADAM17, ADAM18, ADAMTS5, XPNPEP1, CAD, YME1L1, AMZ2, PSMD14, and COPS6. (Figure 2A) Afterward, functional enrichment, including GO and KEGG analysis, was performed. The results of GO analysis indicated that the prognostic MMPs-related DEGs were mainly involved in extracellular matrix degradation, collagen catabolic process, membrane protein ectodomain proteolysis, membrane protein proteolysis, and extracellular matrix organisation (Figure 2B). The KEGG analysis result also suggested that the prognostic MMPs-related DEGs were involved in alanine, aspartate, and glutamate metabolism, Notch signaling pathway, PAR signaling pathway, and IL-17 signaling pathway (Figure 2C).

Figure 2. Identification of prognostic MMPs-related DEGs in hepatocellular carcinoma. (A) Prognostic MMPs-related DEGs are identified from MMPs-related genes through univariate Cox regression analysis. (B) GO enrichment analysis of prognostic MMPs-related DEGs, including biological processes (BP), molecular functions (MF), and cellular components (CC). (C) KEGG enrichment analysis of prognostic MMPs-related DEGs. (D) HCC patients were divided into two clusters according to the consensus clustering matrix (k = 2).

Consensus clustering analysis was then performed based on the expression levels of prognostic MMPs-related DEGs. HCC patients could be divided into the highest intra-group correlation and the lowest inter-group correlation (Figure 2D).

The OS and FPS rates of the two clusters were evaluated, revealing that cluster 1 demonstrated a significantly superior survival status in comparison to cluster 2 (P<0.001) (Figure 3A). The levels of prognostic MMPs-related DEGs were also evaluated in both clusters, and cluster 1 showed lower expression of these genes compared to cluster 2 (P=0.017) (Figure 3B). CIBERSORT was used to determine if the harmful effects of prognostic MMPs-related DEGs in HCC resulted from the infiltration of immune cells. This analysis revealed 22 different immune cell profiles from Cluster 1 and Cluster 2. According to the CIBERSORT analysis, the number of M0 macrophages in cluster 2 showed a substantial decrease when compared with cluster 1 (Figure 3C).

Figure 3. Clinicopathological characteristics of the two clusters classified by prognostic MMPs-related DEGs. (A) Kaplan–Meier curves for the OS and PFS in the two clusters. (B) Comparison of the expression of survival-related MMPs genes in the two clusters. (C) CIBERSORT analysis of infiltrating immune cells in the two clusters.

Construction and validation of the MMPs-related prognostic signature

Subsequently, a LASSO Cox regression analysis was carried out on prognostic MMPs-related DEGs, using the one standard error (SE) approach and 10-fold cross-validation. Consequently, nine prognostic MMPs-related DEGs were determined to be characteristic risk genes. Ultimately, the MMPs-related prognostic signature was formulated using HCC samples from the TCGA-LIHC cohort. The risk score was determined by linearly combining gene expression levels with their respective regression coefficients. The risk score for the correlation coefficient between MMPs-related prognostic signature can be calculated using the following formula: risk score = (0.0086×expression level of RNPEPL1) + (0.0009× expression level of ADAM15) + (0.5255×expression level of ADAM18) + (0.110× expression level of ADAMTS5) + (0.0126× expression level of CAD) + (0.0018×expression level of YME1L1) + (0.0025×expression level of AMZ2) + (0.0114×expression level of PSMD14) + (0.0011×expression level of COPS6) (Figure 4A, 4B).

Figure 4. Construction and validation of the MMPs-related prognostic signature. (A) LASSO regression with tenfold cross-validation found nine prognostic genes using the minimum λ. (B) LASSO coefficient profiles of nine prognostic genes of HCC. (C, D) Heatmap, PCA plot, distribution, survival status, Kaplan-Meier curves for OS and ROC curves demonstrated the predictive efficiency of the risk score in entire, training, and testing cohorts. (E) Validation of the MMPs-related prognostic signature in the training cohort of TCGA-LIHC. (F) The validation of the MMPs-related prognostic signature in the ICGC-LIRI-JP cohort. (G) Univariate analysis of risk score and clinicopathological characteristics in the TCGA-LIHC training and testing cohorts. (H) Multivariate analysis of risk score and clinicopathological characteristics in the TCGA-LIHC training and testing cohorts.

In the entire, training and testing cohorts of the TCGA-LIHC, the heatmap displays the variation in the expression of genes between the high and low-risk groups. According to the median risk score, HCC patients were divided into high-risk and low-risk groups. The TCGA-LIHC cohort showed that the low-risk group had a considerably extended OS duration compared to the high-risk group. The results of PCA indicated a remarkable ability of the risk genes to discriminate between the two groups. The distribution of risk scores indicated a significant increase in the number of fatalities in the high-risk group as compared to the low-risk group. The results of the TCGA-LIHC cohort and the TCGA-LIHC training and testing cohorts showed that patients in the low-risk group had significantly longer OS than patients in the high-risk group (P<0.001). (Figure 4C, 4E) In addition, additional verification was performed in the ICGC-LIRI-JP cohort, and it was observed that the low-risk group had significantly longer OS than the high-risk group (P<0.001). To examine the reliability of the model, the AUC values for OS, which is a time-dependent measure, were calculated by determining the area under the ROC curves. The respective AUC values for the 1-, 3-, and 5-year OS were 0.677, 0.716, and 0.669 (Figure 4F).

Independent predictive value of the risk model

To evaluate the individual impact of the model on HCC prognosis, univariate and multivariate Cox regression analyses were performed in both TCGA-LIHC training and testing cohorts. This analysis involved considering clinical features and the risk score. For both the training and testing cohorts, univariate Cox regression analysis demonstrated that a high-risk score was an independent predictor of adverse patient survival (P<0.001, HR=3.053, 95% CI: 2.071-4.502 and P<0.001, HR=3.316, 95% CI: 1.977-5.563, respectively (Figure 4G).

The outcome was anticipated and validated through the use of multivariate Cox regression analysis, which additionally demonstrated the use of the risk model to be an independent prognostic factor for HCC patients in both the training cohort (P<0.001, HR=2.690, 95% CI: 1.766-4.097) and the testing cohort (P=0.002, HR=2.467, 95% CI: 1.409-4.317), irrespective of other clinical factors (Figure 4H).

To further demonstrate the importance of the prognostic signature associated with MMPs, the ROC results indicate that risk scores can achieve a tremendous overall advantage compared to clinicopathological features in both the TCGA-LIHC and ICGC-LIRI-JP cohorts. The analysis involved examining the correlation between clinicopathological features and risk scores (Figure 5A). In the TCGA group, the risk scores were found to be associated with the T stage and tumor stage (P<0.001), whereas no association was found with age, sex, and grade (P>0.05) (Figure 5B). Furthermore, the associated risk score increases as the tumor stage rises (Figure 5C). The heatmap displays the variation in the expression of the prognostic signature related to MMPs across different clinical characteristics and risk categories within the TCGA-LIHC cohort (Figure 5D). Combining five prognostic factors generated a nomogram to predict the OS at 1-, 3-, and 5-year using the TCGA-LIHC data (Figure 5E). The calibration curves for the prediction of OS at 1-, 3-, and 5-year are in good agreement with the observed values (Figure 5F).

Figure 5. Relationship between risk score and clinicopathological features. (A) Multi-index ROC curve of the MMPs-related prognostic signature in TCGA-LIHC and ICGC-LIRI-JP cohorts. (B) Comparison of the predictive role of risk score for multiple clinicopathological features. (C) Comparison of risk scores for multiple clinicopathological features. (D) Summarized heatmap of the distribution of clinical characteristics and the MMPs-related prognostic signature in TCGA. (E) Nomogram for the quantitative prediction of 1-, 3-, and 5-year survival. (F) Calibration plots for predicting 1-, 3-, and 5-year survival. (G) The results of GO annotation for low-risk and high-risk groups by GSEA. (H) The KEGG annotation results of low-risk and high-risk groups by GSEA.

Functional analysis of MMPs-related prognostic signature

GSEA was performed to identify the biological processes and signaling pathways in which the MMPs-related prognostic signature was enriched. The results showed that the MMPs-related prognostic signature in the high-risk group was mainly enriched in embryonic morphogenesis, external encapsulating structure organization, morphogenesis of an epithelium, ossification, and skeletal system development. In contrast, the low-risk group was significantly enriched in xenobiotic catabolic process, high-density lipoprotein particle, arachidonic acid monooxygenase activity, aromatase activity, and oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen reduced flavin or flavoprotein as one donor and incorporation of one atom of oxygen (Figure 5G). Results from GSEA also indicated that the MMPs-related prognostic signature was mainly enriched in pathways including ECM receptor interaction, cytokine-cytokine receptor interaction, dilated cardiomyopathy, focal adhesion, and neuroactive ligand-receptor interaction in the high-risk group. On the other hand, the group with low risk primarily had an abundance of metabolic pathways such as primary bile acid biosynthesis, beta-alanine, glycine serine and threonine, linoleic acid, and tryptophan metabolism (Figure 5H).

Analysis of immune infiltration in high- and low-risk groups

According to the analysis of ssGSEA, it was found that the group at high risk showed increased levels of expression for four immune cells, specifically aDCs, iDCs, macrophages, and Treg. The low-risk showed increased levels of three immune cells: mast cells, neutrophils, and NK cells (Figure 6A). Moreover, the high-risk group exhibited enhanced functionality of immune cells, such as APC co-stimulation, CCR, Parainflammation, and MHC class I, compared to the low-risk group (Figure 6B). Using CIBERSORT analysis, 22 immune cell profiles were detected in the low and high-risk groups. The results showed that the high-risk group had a higher percentage of M0 and M2 macrophages, while the low-risk group had larger CD8 T and naive B cells (Figure 6C, 6D). Spearman correlation analysis revealed a significant positive correlation between 9 MMPs-related prognostic genes and CD4 memory-activated T cells and M0 macrophages (Figure 6E).

Figure 6. Comparison of immune infiltration and drug sensitivity analysis. (A) The scores of 16 immune cells were detected by ssGSEA analysis. (B) The scores of 13 immune-related functions were detected by ssGSEA analysis. (C, D) 22 types of immune cells were identified by CIBERSORT analysis. (E) Correlation analysis between 22 types of immune cell proportions and MMPs-related prognostic signature. (F) The risk score of various immune infiltration subtypes. (G) Sensitivity of different chemotherapeutic agents in different risk groups.

Additionally, we examined the association between immune infiltrations and risk score to elucidate the role of the MMPs-related prognostic signature in the immune microenvironment. Only C1 to C4 subtypes were detected in patients with HCC, as C5 and C6 subtypes were absent. The results of our study suggest a strong association between immune-infiltrating subcategories, specifically the C1 and C2 subcategories, and elevated risk scores in the TCGA-LIHC cohort. This implies that the MMPs-related prognostic signature might influence the presence of immune infiltrates in individuals diagnosed with HCC (Figure 6F).

Analysis of drug sensitivity

Afterward, we examined the responsiveness of various chemotherapy drugs in distinct risk categories. We observed that the IC₅₀ levels of sorafenib, cytarabine, cisplatin, dihydro rotenone, fludarabine, gemcitabine, irinotecan, oxaliplatin, ribociclib, and mitoxantrone exhibited greater values within the high-risk group. While in the low-risk group, the IC₅₀ values of osimertinib, lapatinib, gefitinib, dasatinib, and 5-fluorouracil were elevated (Figure 6G and Supplementary Figure 1).

Validation of single-cell RNA-seq

Seurat was employed for dimensionality reduction (Supplementary Figure 2) and clustering of the single-cell data, while SingleR was employed for cell-type annotation (Figure 7A, 7B). In the single-cell dataset, the expression of 8 genes belonging to the MMPs-related prognostic signature was identified. These genes include YME1L1, AMZ2, RNPEPL1, ADAM15, PSMD14, COPS6, ADAMTS5, and CAD. Among these genes, YME1L1 and ADAM15 were predominantly expressed in clusters of endothelial cells, while AMZ2, PSMD14, and COPS6 were predominantly expressed in clusters of hepatocytes. Additionally, YME1L1 and COPS6 were the main genes expressed in clusters of monocytes (Figure 7C, 7D).

Figure 7. Single-cell RNA-seq profiling in hepatocellular carcinoma. (A) t-SNE and UMAP plots colored by different cell clusters. (B) The cell types are identified by marker genes. (C) Expression of MMPs-related prognostic signature in each cluster.