Multi-omics analysis of pyroptosis regulation patterns and characterization of tumor microenvironment in patients with hepatocellular carcinoma

Data collection

Gene expression profile data and clinical information for HCC were collected from TCGA (https://portal.gdc.cancer.gov/). The clinical data included age, sex, and survival status (Table S1). We obtained 369 tumor and 50 normal tissue samples. The copy number variation data for HCC were downloaded simultaneously. In addition, the Chip data set of GSE76426 (Grinchuk et al., 2018) from the GEO database was analyzed as a validation dataset and included 115 primary tumors (PTs) tissue samples and 53 adjacent non-tumor tissues (ANTTs) samples derived from 115 HCC patients, The dataset was obtained using the GPL 10558 Illumina Human HT-12 V4.0 expression beadchip sequencing platform. All the gene sets were derived from Homo sapiens. The overview of the workflow is shown in Fig. 1.

Differentially expressed genes of HCC

To analyze the influence of gene expression values on the occurrence and development of HCC, differential gene analysis was performed on normal and tumor samples from TCGA dataset using the DESeq2 R package (Love, Huber & Anders, 2014). To define the threshold of DEGs, we applied the absolute value of log₂Fold-change > 2 and Padj < 0.05; cut-off values of logFC > 2 and Padj < 0.05 and logFC < −2 and Padj < 0.05 were used for upregulated and downregulated genes, respectively. To uncover the influence of pyroptosis-related genes on biological functions, we extracted data from the GOBP_PYROPTOSIS dataset of the Genecards database (http://www.genecards.org) (Safran et al., 2010) and REACTOME_PYROPTOSIS dataset of the MSigDB (Liberzon et al., 2015) database. We then examined the intersection of pyroptosis-related genes and DEGs in HCC; the results were represented in the form of Venn diagrams. A volcano plot was constructed to show HCC-related and pyroptosis-related DEGs.

Construction of risk model for HCC

To analyze the effect of pyroptosis-related DEGs on the prognosis of patients with HCC, the expression and survival data of patients from TCGA database were analyzed. Prognostic risk genes for HCC were evaluated using univariate Cox proportional hazards regression analysis. After the prognostic risk genes for HCC were included in the model, we used the least absolute shrinkage and selection operator (LASSO) algorithm for dimension reduction and obtained prognostic feature genes. Each normalized gene expression value was weighted by penalty coefficients using LASSO-Cox analysis, and a risk score formula was established. Patients were then divided into high- and low-risk groups based on the mean risk score calculated as follows:

$$riskScore={\displaystyle \sum _{i}C}oefficient(risk{\text{gene}}_{\text{i}})*\text{mRNA\hspace{0.17em}Expression\hspace{0.17em}}(\text{risk\hspace{0.17em}}{\text{gene}}_{\text{i}}).$$

To verify the accuracy of the risk score, we used the GSE76427 dataset, which included patient gene expression and clinical data, to group patients and conduct a statistical analysis of the differences in survival.

Clinical prediction model based on risk model

We used a combination of the risk score combined and clinicopathological characteristics to assess the personalized prognosis of patients. Univariate and multivariate Cox regression analyses of risk score and clinicopathological characteristics were used to assess the predictive power of OS. Finally, a clinical predictive nomogram model was constructed by incorporating the risk score model and important clinicopathological parameters. To quantify the discriminative performance, we compared the nomogram-predicted probability with the observed actual survival rate and used a calibration curve to evaluate its nomogram performance.

Differentially expressed genes between the high-risk and low-risk groups

To analyze the effect of the risk score on HCC, DEG analysis of the samples that were divided into high- and low-risk groups from TCGA dataset was performed using the R package DESeq2 (Love, Huber & Anders, 2014). Subsequently, the significant differential genes were screened. The threshold values to determine the DEGs were as follows: the absolute values of the log₂fold change (log₂FC) > 2 and Padj < 0.05 were defined as upregulated DEGs and log₂FC < −2 and Padj < 0.05 were defined as downregulated DEGs. The results were then represented as a volcano plot.

We obtained the dataset of HCC cell line-drug effects from the GDSC database (https://www.cancerrxgene.org/) (Yang et al., 2013). The gene expression data of high- and low-risk patients from TCGA database were subjected to drug sensitivity analysis using the R package OncoPredict (Maeser, Gruener & Huang, 2021). Subsequently, the sensitivity differences of the therapeutic drugs for HCC were compared between the high- and low-risk patients. In addition, to determine if the genes in the high- and low-risk patients with HCC showed copy number variation, GISTIC 2.0 in Genepattern (https://cloud.genepattern.org/) (Reich et al., 2006) was used to analyze the copy number variation of patients with HCC from TCGA database.

Assessment of biological characteristics

We used GO, KEGG (Ashburner et al., 2000; Kanehisa & Goto, 2000), GSEA, and GSVA (Hänzelmann, Castelo & Guinney, 2013) to analyze the biological characteristics of patients of different risk groups. GO annotation analysis and KEGG pathway enrichment analysis were performed on the DEGs using the clusterProfiler package (FDR < 0.05) (Yu et al., 2012). To investigate whether there were biological process differences between these two patient groups, we downloaded “c5.go.v7.4.entrez.gmt” and “c2.cp.kegg.v7.4.entrez.gmt” as reference gene sets from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb) based on the gene expression profile data sets of patients with HCC. GSEA was performed, and the results were visualized using the clusterProfiler R package (P < 0.05). To study the variation in the biological processes of high-risk samples compared to that in those of low-risk samples, datasets of gene expression profiles based on TCGA were used to perform GSVA using the GSVA package in R. The reference gene set (“h.all.v7.4. symbols.gmt”) was downloaded from MSigDB. We then calculated an enrichment score for each sample in each pathway and examined the correlation between the enrichment score for patients and the risk score (P < 0.05).

Identification and correlation analysis of tumor-infiltrating immune cells

Analysis of immune cell infiltrates in tissues is a key component in the characterization of diseases and in the prediction of disease prognosis. The stomal and immune cell contents of high- and low-risk samples from TCGA dataset were estimated using the R package “estimate” (Yoshihara et al., 2013). The correlations between the risk score and ESTIMATE score were calculated. We evaluated the proportion of 22 immune cells in the immune microenvironment of high- and low-risk samples using the CIBERSORT algorithm (Newman et al., 2019). We set the number of permutations to 1,000 and used the Wilcoxon test to calculate the differences in proportions of immune cells between the high- and low-risk groups (P < 0.05).

Construction of protein-protein interaction (PPI) network

We used the STRING (Szklarczyk et al., 2019) dataset to construct PPIs of differentially expressed pyroptosis-related genes with a combined score of >400. Cytoscape (v3.7.0) was used for visualization of the PPI network, and Clue GO (Bindea et al., 2009) was used for functional annotation. Next, we used the cytoHubba plug-in to obtain hub genes (Chin et al., 2014) using the MCC method in the PPI network.

RNA extraction and real-time quantitative PCR assay

Total RNA was extracted from liver cancer and adjacent tissues (samples of liver cancer and adjacent tissues were collected from volunteers who underwent liver cancer surgery at the Second Affiliated Hospital of Dalian Medical University between January and December 2012, and each subject signed an informed consent form. All patient samples were approved by the Ethics Committee of the Second Affiliated Hospital of Dalian Medical University. No. 134.) using the M5 HiPer Universal RNA Mini Kit (Mei5 Biotechnology Co. Ltd, Beijing, China). We then used the PerfectStart® Uni RT&qPCR Kit (TransGen Biotech, Beijing, China) for reverse transcription into cDNA and performed qRT-PCR. The primers were as follows: (MDK) forward: 5′-GCAACTGGAAGAAGGAGT-3′, and MDK reverse: 5′-TGGCACTGAGCATTGTAG-3′.

Immunohistochemistry

Immunohistochemistry was performed on liver cancer and adjacent tissues (the source of sample is the same as that of qRT-PCR experiment). We performed dewaxing with xylene, followed by a gradient dehydration with alcohol, and then inactivated the endogenous peroxidase with hydrogen peroxide. The sections were then blocked and treated with primary antibodies against MDK (anti-Midkine; Abcam, Cambridge, UK). The antibody was diluted with PBS at a ratio of 1:600, and then incubated for 12 h at 4 °C. To each film, 100 µL of working fluid (PV-9000; Beijing Zhongshan Golgen Bridge Biological Technology Co. Ltd, Beijing, China) was added and then incubated at 37 °C for 20 min. The solution was dehydrated again after dyeing, and the film was sealed.

Statistic analysis

All data were processed and analyzed using R software (version 4.1.1). Continuous variables were compared between the two groups. For normally distributed variables, we used independent t-tests; for non-normally distributed variables, we used the Mann-Whitney rank-sum test (Li et al., 2022). Correlations between different genes were analyzed using the Pearson correlation coefficient calculation. Survival analysis was performed using the R survival package. Survival differences are shown using Kaplan–Meier survival curves. The significance of the difference in survival between the two groups was assessed using the log-rank test. Univariate and multivariate Cox regression analyses were used to identify independent prognostic factors. All statistical P-values were two-sided, and P < 0.05 was considered statistically significant.

Differentially expressed pyroptosis-regulated factors and the effect on their biological processes

To reveal the biological differences between HCC and normal samples from transcriptomes, we performed differential gene expression analysis in HCC and normal samples. After screening, 8,958 genes were defined as significant DEGs, comprising 7,030 upregulated and 1,928 downregulated genes (Fig. 2A). In total, 197 pyroptosis-related genes were obtained from the GeneCards database. The intersection of DEGs and pyroptosis-related genes revealed 37 differentially expressed pyroptosis-related genes (22 upregulated and 16 downregulated; Figs. 2B, 2C).

Functional enrichment was performed for differentially expressed pyroptosis-related genes (Fig. 2D). The results showed that the differential genes were related to biological processes (Fig. 2E) including digestion, maintenance of gastrointestinal epithelium, digestive system process, epithelial structure maintenance, and O-glycan processing; cellular components (Fig. 2F) including extracellular space, extracellular region, endosome membrane, and plasma membrane; and molecular functions including interleukin-1 receptor binding, lipoteichoic acid binding, G-protein coupled adenosine receptor activity, lipopolysaccharide binding, and identical protein binding (Fig. 2G). Furthermore, KEGG pathway enrichment analysis indicated that the differentially expressed pyroptosis-related genes were significantly enriched in the pertussis, Salmonella infection, toll-like receptor signaling, and tuberculosis pathways (Fig. 2H).

Construction of risk models and prognostic analysis

To analyze the effects of differentially expressed pyroptosis-related genes on the prognosis of patients with HCC, we used univariate Cox regression analysis and identified 45 prognostic risk genes of HCC. The prognosis-associated risk genes of HCC were included in the LASSO-Cox analysis to choose 12 genes with the best prognostic value. We subsequently performed a correlation analysis of the expression of the 12 characteristic genes, and the results showed a high level of correlation between these genes (Fig. 3A). The penalty coefficients of characteristic genes were calculated based on LASSO-Cox analysis (ADRA1A: −0.023, GPRC6A: 0.037, MAGEA8: 0.034, MDK: 0.024, MMP9: 0.019, NETO2: 0.027, OLFM3: 0.031, PYCR1: 0.020, RET: 0.029, and TNFRSF4: 0.065). The risk score was calculated by multiplying the expression level of a gene with its corresponding coefficient and then adding them together. Patients with HCC were then divided into high- and low-risk groups based on their mean risk scores. In addition, principal component analysis revealed that the 12 characteristic genes could better distinguish patients in both TCGA and GSE76427 datasets (Figs. 3B, 3C). Kaplan–Meier analysis revealed worse OS scores in patients with high risk scores (log-rank P < 0.0001, Fig. 3D). The risk score of each patient with HCC in GSE76427 was then computed, and the Kaplan–Meier analysis revealed that patients with a high-risk score showed worse OS (log-rank P < 0.05, Fig. 3E).

Construction of clinical prediction model based on risk score

We evaluated the prognostic impact of the risk score in patients with HCC. Univariate and multivariate Cox analyses revealed that the risk score was an independent prognostic factor for patients with HCC (Figs. 4A, 4B). Different clinicopathological features were combined into a risk score to generate a prediction model nomogram for predicting the OS of patients with HCC (Fig. 4C). The calibration curve showed a good prediction of 1-, 2-, and 3-year survival for patients using the model (Figs. 4D–4F). Concurrently, the results of time-dependent receiver operating characteristic curve analysis showed that the predicted percentages of 1-, 3-, and 5-year survival were 77.8%, 79.8%, and 82.4%, respectively (Fig. 4G). The decision curve analysis supported the model, which could offer greater prognostic gains for patients (Fig. 4H).

Divergence analysis between high-risk group and low-risk group

Owing to the clear distinction in survival between the high- and low-risk groups, we performed differential expression analysis of the expressed genes in the high- and low-risk groups. In total, 349 DEGs were identified. Of these, 51 genes were upregulated and 298 were downregulated (Fig. 5A). Concurrently, the DEGs were subjected to dimensionality reduction using principal component analysis. The results showed a significant difference between the high- and low-risk groups (Fig. 5B). We then analyzed the effect of differentially expressed mRNA on the biological functions associated with the high- and low-risk groups. GO biological function annotation analysis of the DEGs showed that these genes were mainly enriched in digestion, maintenance of gastrointestinal epithelium, digestive system process, epithelial structure maintenance, O-glycan processing, and other biological processes (Fig. 5C); catenin complex, GABAergic synapse, Golgi lumen, dendrite membrane, apical part of cell, and other cellular components (Fig. 5D); and ligand-gated anion channel activity, chloride channel activity, hormone activity, anion channel activity, and other molecular functions (Fig. 5E; Table S2). Additionally, KEGG pathway enrichment analysis of the DEGs showed that these genes were mainly enriched in pathways associated with neuroactive ligand-receptor interaction, taste transduction, pancreatic secretion, taurine and hypotaurine metabolism, and the cAMP signaling pathway (Fig. 5F; Table S3). We also analyzed the copy number differences between the high- and low-risk groups. The results revealed that the copy number amplification of patients assigned to the high-risk group was significantly lower than that of those assigned to the low-risk group, and the high-risk group included more chromosomal deletions (Fig. 6).

To analyze the influence of risk score on the treatment of HCC, we integrated the gene expression data of high- and low-risk patients in TCGA dataset with drug sensitivity for HCC from the GDSC database. We identified drugs that were more sensitive in high- and low-risk patients and had lower IC₅₀ values. Among the drugs currently used for the treatment of HCC, IGFR_3801, tanespimycin, docetax, and lestaurtinib were more sensitive in high-risk patients (Fig. 7A), whereas belinostat, ispinesib mesylate, shikonin, and midostaurin were more sensitive in low-risk patients (Fig. 7B). Among the currently used drugs for other cancers, trametinib, mitomycin-C, rapamycin and obatoclax mesylate were more sensitive in high-risk patients (Fig. 7C), whereas pevonedistat_1529, methotrexate_1008, foretinib_308, and AZD8055_1059 were more sensitive in low-risk patients (Fig. 7D).

Next, we performed GSEA of the gene expression data of the high- and low-risk groups. We found that the patients in the high- and low-risk groups showed significant differences mainly in the following biological processes: inhibitory extracellular ligand-gated ion channel activity, limb bud formation, condensed chromosome outer kinetochore, spindle elongation, regulation of chloride transport, and other processes were inhibited, whereas very low-density lipoprotein particle remodeling, triglyceride-rich lipoprotein particle remodeling, complement activation lectin pathway, platelet dense granule lumen, alcohol dehydrogenase were activated (Figs. 8A and 8B). Simultaneously, we found that primary bile acid biosynthesis, fatty acid metabolism, complement and coagulation cascades, glycine serine and threonine metabolism, renin angiotensin system, and other pathways were activated, whereas some pathways, such as cell cycle, ECM receptor interaction, neuroactive ligand receptor interaction, gap junction, and axon guidance, remained inactivated (Figs. 8C and 8D; Tables S4 and S5).

We further used GSVA to calculate the DEGs of high- and low-risk patients from TCGA database to determine their effect on biological characteristics and oncogenic signaling pathways. The results revealed that the risk score was closely correlated with the activity of multiple oncogenic signaling pathways (Fig. 8E). Among them, there was a significant negative correlation between risk score and hallmark KRAS signaling DN, hallmark pancreas beta cells, hallmark coagulation, hallmark xenobiotic metabolism, and hallmark bile acid metabolism, whereas there was a significant positive correlation between risk score and hallmark G2M checkpoint, hallmark E2F targets, hallmark MYC targets V2, hallmark DNA repair, and hallmark MYC targets V1.

Immune characterize differences in high-risk and low-risk groups

We analyzed the differences in immune cell content in the high- and low-risk groups (Fig. 9A), and the high- and low-risk groups showed significant differences in multiple immune cells (Fig. 9B). The risk score was positively correlated with B cells memory, dendritic cells resting, macrophages M0, and regulatory T cells (Tregs), and negatively correlated with resting mast cells, monocytes, and T cells gamma delta (P < 0.05, Fig. 9C). Additionally, we calculated the correlation between immune cells and patients in the high- and low-risk groups. We found that there was a positive correlation between T cells CD4 memory resting/dendritic cells resting in patients in the high-risk group (Fig. 9D), while there was a negative correlation between B cells memory, dendritic cells resting, neutrophils, Tregs, macrophages M0 in patients in the low-risk group (Fig. 9E).

We also analyzed the correlations between the risk score of HCC and the immune score, stromal score, tumor purity, and ESTIMATE score, among which there was a significant positive correlation between the risk score and stromal score (P < 0.05, Fig. 10A). We constructed a PPI network of DEGs in patients in the high- and low-risk groups (Fig. 10B). Visualization of the PPI network was carried out with Cytoscape, which included 159 DEGs and 272 PPI pairs in total. The top five genes that interacted closely with other DEGs were MUC5AC (interacted with 16 DEGs), MUC5B (interacted with 16 DEGs), MUC1 (interacted with 15 DEGs), CFTR (interacted with 14 DEGs), and CHGA (interacted with 13 DEGs). In addition, DEGs related to the PPI network were mainly enriched in functions related to calcium-ion regulated exocytosis, digestive system processes, and cochlear development (Fig. 10C). The hub genes of the PPI network were screened using the cytoHubba plugin, and these genes are more likely to be the key genes causing differences between the high- and low-risk groups (Fig. 10D).

High expression of MDK in hepatocellular carcinoma

From the RT-qPCR results (Fig. 11), it was clear that the expression of MDK in HCC is much higher than that in adjacent tissues. The immunohistochemical analysis (Fig. 12) also showed that MDK was highly expressed in HCC tissues compared to that in normal liver tissues.