Analysis of clinical significance and molecular characteristics of methionine metabolism and macrophage-related patterns in hepatocellular carcinoma based on machine learning

2.1 Data acquisition and processing

mRNA expression and patient clinical data were downloaded from the gene expression omnibus (GEO) database (GSE14520, GSE54236, GSE76427, and GSE116174), The Cancer Genome Atlas (TCGA) database [15] (TCGA-LIHC dataset), Genotype-Tissue Expression (GTEX) database (GTEX liver dataset) [16], and International Cancer Genome Consortium (ICGC) database (ICGC-JP dataset) [17]. Batch calibration of different datasets was accomplished using the “SVA” (3.1.1) R package. Gene difference analysis was performed in the form of count data and corrected using the limma software package. Downstream function analysis utilized the form of TPM data. Differential expression analysis was performed using the “LIMMA” (3.54.1) R software package. The protein expression data of HCC patients were downloaded from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database [18].

2.2 Identification of methionine metabolism- and macrophage-related gene modules

The methionine metabolic pathway score was based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) downloaded dataset (Table S1) and the “GSVA” (1.46.0) R package for quantitative scoring [20]. Tumour-associated macrophages in the HCC group were quantified using the previously published “EPIC” and “CIBERSORT” algorithms [20]. The “WGCNA” (1.68) R package was used for weighted correlation network analysis (WGCNA), and an appropriate soft threshold was chosen to cluster genes with similar coexpression in the same module [19]. Methionine metabolism and macrophage infiltration data were combined with the above modules to identify related gene clusters.

2.3 Recognition of methionine metabolism- and macrophage-related gene expression patterns in hepatocellular carcinoma

The aforementioned identified gene sets were subjected to prognostic correlation analysis, and the critical value of p was less than 0.05. Then, the K-means algorithm was used to identify the above gene expression patterns in GEO meta-datasets using the “ConsensusClusterPlus” (1.58.0) R software package [19]. The same algorithm was used to verify the above results in the TCGA-LIHC dataset.

2.4 Functional prediction analysis of the different gene patterns

In the different gene patterns identified above, the differences in signalling pathways and immune cell infiltration among different subgroups were evaluated by using the “clusterProfiler” (4.2.2) R package and “CIBERSORT” algorithm 20. The “survival” (3.5-0) R package was used for survival analysis [19].

2.5 Construction of the diagnostic and prognostic model

The aforementioned 47 related genes used for genotyping were further analysed for differences, and the cut-off value was an absolute value of logFC greater than 1. Based on 5x cross validation, a diagnosis model was constructed by using random forest and logistic regression, and a ROC curve was used for diagnosis efficiency. We verified the above built models in the independent TCGC-LIHC, ICGC-JP and CPTAC datasets. The prognostic model was used to analyse the prognosis of the aforementioned 47 genes in the TCGA-LIHC dataset, and a total of 21 prognostic genes were identified. Screening of important variables was accomplished with lasso regression analysis. Multivariate Cox regression analysis was used to construct a prognostic model, and the aforementioned model was validated in the TCGA-LIHC dataset and the independent ICGC-JP dataset [21].

2.6 Evaluation of drug effects

The CTRP2.0 database, which contains the sensitivity data for 481 compounds over 835 CCLs, and the PRISM database, which contains the sensitivity data for 1448 compounds over 482 CCLs, were used to evaluate the efficacy of compounds according to K-nearest neighbour (KNN) imputation. The top and bottom 20% of the high and low subgroups were used as cut-off points. The specific methods were described by Shixue Dai et al. [22].

2.7 Immunohistofluorescence

Samples were prepared through gradient dehydration of paraffin sections and EDTA antigen thermal repair. Primary antibody (MAT1A, affinity cat. DF12210, dilution 1:200 with TBST; and CD68, Proteintech cat. 3A9A7, dilution 1:300 with TBST) was incubated overnight. After three washes with PBS, immunofluorescence secondary antibody (FITC Anti-Mouse, cat. Ab6785, dilution 1:500; Alexa Fluor ® 647 Anti-Rabbit, cat. ab150083, dilution 1:500) was incubated at room temperature for 1Â h. The nuclei were stained with DAPI. All tissues were observed under a fluorescence microscope with the same intensity light source and exposure time for quantitative comparison.

2.8 Statistical analysis

SPSS (22.0), R (4.1.0) and Python (3.4) software were used for statistical analysis. Unless otherwise indicated, a P value 0.05 was considered statistically significant. Two-tailed Student’s t tests were used to test the differences between different groups.

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

3.1 Screening methionine- and macrophage-related gene modules in hepatocellular carcinoma

Based on the four large human HCC datasets (GSE14520, GSE54236, GSE76427, and GSE116174), which were named the GEO meta-dataset in the GEO database, we used GSVA to score the methionine metabolic pathway of each HCC sample as previously reported and used EPIC and CIBERSORT, which are currently the most widely used, to quantify tumour tissue macrophages (Fig. 1A). Furthermore, we used the GEO meta-dataset to identify differentially expressed genes. The critical values were an absolute value of logFC greater than 0.585 and P less than 0.05. A total of 1314 differentially expressed genes were identified, including 650 downregulated genes and 664 upregulated genes (Fig. 1B and Table S2). Based on WGCNA, we divided the above 1314 genes into six modules and conducted coexpression analysis with macrophages and methionine metabolic activity. It was found that the turquoise module was the most significantly related to the above, so we focused on it in subsequent research (Fig. 1C).

3.2 Identification of two HCC patterns and correlation with prognosis

In the survival analysis, with a cut-off value of p less than 0.05, we identified a total of 47 prognosis-related genes in the aforementioned gene module for downstream analysis (Table S3). Based on the GEO meta-dataset, we used the k-means algorithm to cluster and found that HCC can be divided into two patterns based on macrophage- and methionine metabolism-related genes (Fig. 2A). PCA based on the expression of all genes also showed that the above clustering was obvious (Fig. 2B). Survival analysis showed that the prognosis of Cluster 1 was significantly better than that of Cluster 2 (Fig. 2C). In the independent TCGA-LIHC dataset of HCC, we verified the above results, which suggests the reliability of our above typing (Fig. 2D–F).

3.3 Molecular characteristics of the two tumour patterns

The above patterns were matched with clinical traits for correlation analysis. We found that pattern 2, with a worse prognosis, had a later TNM and BCLC stage and a larger tumour volume. These results internally supported the aforementioned survival analysis results (Fig. 3A). Interestingly, in pattern 2, there was less tumour-associated macrophage infiltration and weaker methionine metabolic activity (Fig. 3B and C). KEGG signalling pathway analysis showed that pattern 2 had weaker fatty acid metabolism and PPAR pathway activity and stronger cell cycle activity (Fig. 3D and Table S4). Similar results were obtained based on TCGA dataset analysis (Fig. 3E–G). Using immunofluorescence double labelling, our own HCC cohort also supported our clustering results (Fig. 3H).

3.4 Construction of the HCC diagnosis model based on the above methionine- and macrophage-related genes

Although the pathological examination of HCC is the gold standard for the diagnosis of HCC, it is limited by the professional skills of pathologists, especially in some smaller medical institutions. We tried to build a diagnostic model of HCC by using the previously identified genes. First, we raised the differential gene expression threshold to the absolute value of logFC = 1, identified 22 candidate genes, and constructed a diagnosis model by using random forest and logistic regression models (Fig. 4A and B). In the GEO meta-dataset of the training set, the random forest and logistic regression models showed similar excellent discrimination efficiency. However, in the test set (TCGA, ICGC and CPTAC datasets), the logistic regression model was significantly better than the random forest model (Fig. 4C–F).

3.5 Construction of a prognostic model for hepatocellular carcinoma based on the above methionine- and macrophage-related genes

We performed univariate Cox regression analysis in the TCGA dataset based on the previously identified methionine- and macrophage-related genes. Twenty-one prognosis-related genes were identified (Fig. 5A). Using lasso regression analysis, we further identified the key variables, and a total of 8 genes were screened and used to construct the prognostic model (Fig. 5B). Using multivariate regression analysis, we constructed a prognostic model for HCC. In the training set GEO meta-dataset and the independent ICGC and TCGA datasets, the ROC curve showed good discrimination efficiency of the model (Fig. 5C–E). Clinical correlation analysis showed that patients in the high-risk group had later BCLC and TNM stages, further supporting the reliability of our model (Table 1).

To make better use of the above prognostic model in the clinic, we conducted univariate and multivariate Cox regression analyses and found that the above risk score and BCLC and TCGA stages were independent prognostic factors (Fig. 6A and B). A nomogram was constructed based on the TCGA staging, BCLC staging and risk score (Fig. 6C). In the GEO meta-dataset, TCGA dataset, and ICGC dataset, the nomogram above showed a better ability to distinguish prognosis (Fig. 6D–F).

3.7 Bioinformatics prediction of drug candidates

To identify candidate drugs for the high-risk group, we comprehensively analysed the PRISM and CTRP2.0 databases and found that four compounds may be effective in the low-risk group based on the KNN algorithm. The results based on PRISM found that six drugs had better responses in the low-risk group (Fig. 7A), and the results based on CTRP2.0 found that three drugs had better responses in the low-risk group (Fig. 7B). These results can provide some guidance for clinical treatment.