Biomarkers linked to M2-like tumor-supporting immune cells for predicting outcomes, understanding immune environment, and finding drugs in liver cancer

The GSE125449↗ single-cell transcriptome profiles of liver cancer was downloaded from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/↗). We selected seven HCC samples from Set 1 for analyses. The “Seurat” package (13) was used for processing scRNA-seq data, including data filtering (cells and genes), normalization, principal component analysis (PCA), and t-distributed stochastic neighbor embedding (t-SNE). The quality control standards referred to the uploader (12). Cell samples with >20% mitochondrial gene expression were filtered. Cells with >700 detected genes and genes detected in >3 cells were reserved. The “DoubletFinder” package (14) was used to remove samples with a doublet rate >0.4%. After cell filtering, the scRNA-seq data of high-quality cells were normalized to find highly variable genes for downstream analyses. Then, PCA was done on highly variable genes to identify significant principal components (PCs). Cell clustering was undertaken on the top-20 PCs using the t-SNE algorithm. The “FindAllMarkers” function was applied to detect the marker genes of each cell cluster. Next, annotation of cell type in different cell clusters was done with the “SingleR” package (15). HCC-related clinical information and gene-expression data were downloaded from The Cancer Genome Atlas (TCGA) database (www.cancer.gov/↗), GEO database (GSE76427↗), and the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/↗), and only HCC samples with complete survival information were retained. Then, the TCGA-LIHC dataset (which contains the survival data and clinical information for 368 HCC patients) was used as the training set. Gene-expression data from TCGA-LIHC were downloaded in the format of fragments per kilobase million and analyzed. Data on progression-free survival (PFS), disease-specific survival (DSS), and disease-free survival (DFS) for TCGA-LIHC were downloaded from UCSC Xena (https://xena.ucsc.edu/↗) (16). The ICGC-LIRI-JP dataset and GES76427↗ dataset contain the survival data, clinical information, and gene-expression data for 232 and 115 HCC patients, respectively. The mRNA-seq data in the ICGC-LIRI-JP dataset were transformed by log₂(x+1), and data in the GSE76427↗ dataset were normalized using the “limma” package (17). Then, the batch effect between the ICGC-LIRI-JP dataset and GSE76427↗ dataset was eliminated using the “sva” package (https://bioconductor.org/packages/sva/↗), so that they were combined into a merged dataset to serve as the test set. A summary of the clinicopathological characteristics of patients in all datasets is shown in Supplementary Table S1.

The relative content of M1 and M2 macrophages in each TCGA-LIHC sample was calculated on CIBERSORTx (https://cibersortx.stanford.edu/↗) (18) using the default signature matrix. The “surv_cutpoint” function of the “survminer” package (https://rdocumentation.org/packages/survminer/↗) was used to calculate the optimal cutoff value to distinguish high- and low-content groups of M1 or M2 macrophages in TCGA-LIHC samples. Survival analyses were carried out using the “survival” package (https://cran.r-project.org/web/packages/survival/index.html↗). Survival between low- and high-M1 (or M2) macrophage-content groups were analyzed and compared by Kaplan–Meier method to ascertain if M1 and/or M2 macrophage content was related to survival from HCC.

After grouping the HCC samples by trait of high or low M2 macrophage content, we analyzed TCGA-LIHC expression data using the Weighted Gene Co-expression Network Analysis (“WGCNA”) package (19) to obtain genes most related to M2 macrophage content. Samples were clustered to ascertain the overall relevance of all samples in the dataset, and outliers were excluded. The soft thresholding power β was chosen based on the lowest power for which the scale-free topology fit index reached a high value. The minimum gene number/module was set to 50 and, finally, 11 modules were generated. Next, we undertook correlation analyses between modules and traits to find the most relevant modules for M2 macrophage content. Finally, the obtained modular genes were intersected with the TAM marker genes acquired from analyses of scRNA-seq data to filter M2-like TAM-related genes.

To obtain M2-like TAM-related genes that could construct a prognostic signature, univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses were carried out. Initially, we wished to uncover the association between signature genes and the prognosis. Hence, after the consensus clustering of HCC samples into different clusters based on expression of signature genes, we analyzed the difference in the prognosis among clusters. And we undertook analyses of the enrichment of function and signaling pathways of signature genes using the Gene Ontology (GO) database (http://geneontology.org/↗) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database (www.genome.jp/kegg/↗), respectively, by employing the “clusterProfiler” package (20). To group the HCC patients, the risk score of each HCC patient in the training set was calculated according to the following formula:

Then, patients were divided into low- and high-risk groups based on the optimal cutoff of the risk score. Kaplan–Meier survival curves and the log-rank test were used to analyze and compare the survival between low- and high-risk groups. Receiver operator characteristic (ROC) curves for 1, 3, and 5 years were plotted using the “survivalROC” package (https://cran.rstudio.com/web/packages/survivalROC/index.html/↗) to evaluate the performance of the prognostic signature. According to the KEGG database, signaling pathways enriched significantly in low- and high-risk groups were analyzed by gene set enrichment analysis (GSEA) software (21). The count of permutations was set to 1000. Significantly enriched pathways were defined as those with P<0.05 and a false discovery rate<0.25. The five most enriched pathways in the high- and low-risk groups, respectively, were obtained. Moreover, the prognostic signature was validated in the test set. Based on the prognostic signature and clinical characteristics of samples, the “rms” package (https://cran.r-project.org/web/packages/rms/index.html↗) was used to construct a nomogram. The performance of the nomogram was evaluated using calibration curves and 1-, 3-, and 5-year ROC curves.

The “GSVA” (22) and “GSEABase” packages (www.bioconductor.org/packages/release/bioc/html/GSEABase.html/↗) were used to analyze differences in scores for immune cells and immune function between high- and low-risk groups. The Tumor Immune Dysfunction and Exclusion (TIDE) score was calculated for each sample in high- and low-risk groups on the TIDE website (http://tide.dfci.harvard.edu/↗) (23). The immunophenoscore of each sample was obtained on The Cancer Immunome Atlas (TCIA) database (https://tcia.at/home/↗) (24).

We wished to further identify new potential targets and more efficacious drugs for HCC treatment. The CellMiner database (https://discover.nci.nih.gov/cellminer/↗) was employed to screen for antitumor drugs whose sensitivity was associated significantly with prognostic genes. The “pRRophitic” package (https://github.com/paulgeeleher/pRRophetic/↗) was used to predict the half-maximal inhibitory concentration (IC₅₀) of different drugs in high- and low-risk groups. The lower the IC₅₀ of a drug, the more efficacious the drug is for treating cancer.

The study protocol was approved by the Human Subjects Committee of Xijing Hospital (Xian, China). All patients provided written informed consent. We collected samples of tumor tissue and adjacent normal tissue from 15 patients with HCC. Detailed clinicopathological information is summarized in Supplementary Table S2. Total RNA from human tissues was extracted using TRIzol^® Reagent (Invitrogen). Then, the RNA was reverse-transcribed into complimentary-DNA using a PrimeScript RT kit (Takara Biotechnology. Shiga, Japan). qPCR was done using SYBR Premix Ex Taq II (Takara Biotechnology) for a real-time PCR detection system (Bio-Rad Laboratories, Hercules, CA, USA). Supplementary Table S3 lists all the primers used in PCR. Expression of genes was normalized to that of glyceraldehyde 3-phosphate dehydrogenase (GAPDH).

Statistical analyses were conducted using R 4.0.3 (R Institute for Statistical Computing, Vienna, Austria). The packages within R used for statistical analyses were as described above. The Kaplan–Meier method was employed for survival analyses. The Wilcoxon test was used to compare differences between two groups. The Kruskal–Wallis test was employed to compare differences among three or more groups. P<0.05 was considered significant.

We wished to further clarify the relationship between macrophages and the HCC prognosis. The CIBERSORTx algorithm was used to calculate the content of M1 and M2 macrophages in TCGA-LIHC samples. Then, HCC patients were divided into high- and low-M1 macrophage-content groups and high- and low-M2 macrophage-content groups. Kaplan–Meier analyses showed no significant difference in survival from HCC between high- and low-M1 macrophage-content groups (Figure 2A), but HCC patients in the low-M2 macrophage-content group had longer survival (Figure 2B), thereby indicating that M2 macrophages had an important role in HCC. Based on this observation, WGCNA was undertaken to identify M2 macrophage-related genes in HCC. First, no outlier was detected in TCGA-HCC (Figure 2C), and 7 was chosen as the optimal soft-threshold power (Figures 2D, E), and 11 modules were identified by WGCNA (Figures 2F, G). Correlation analysis between modules and M2 macrophage content showed the red module to be associated most significantly with high-content M2 macrophages (correlation = 0.32, P<0.001). Thus, 405 genes (Supplementary Table S4) in the red module were selected for downstream analyses.

After quality control, 19,106 genes within 2,719 cells were obtained. The number of genes (nFeature), the sequence count per cell (nCount), and percentage of mitochondrial genes (percent.mt) were displayed in Vlnplots (Figure 3A). Correlation analyses showed that nCount was correlated positively with nFeature (Figure 3B). Then, 2000 variable genes were plotted in a scatter diagram (Figure 3C). Thirty PCs were identified (Figures 3D, E), showing high heterogeneity in HCC cells. The top-20 PCs were selected for t-SNE analyses. According to t-SNE and cell-type annotation, HCC cells were clustered into two groups: 1,226 immune cells and 1,493 non-immune cells (Figure 3F). The immune group was composed of B cells, T cells, and TAMs. The non-immune group included cancer-associated fibroblasts (CAFs), cells with an unknown entity but express hepatic progenitor cell markers (HPC-like cells), malignant cells, tumor-associated endothelial cells (TECs), and unclassified cells (Figure 3G). The 2047 TAM marker genes of immune cells were detected (Supplementary Table S5) and shown in a heatmap (Figure 3H).

After marking the intersection of 2047 TAM marker genes and 405 M2 macrophage modular genes, 127 candidate M2-like TAM-related genes were obtained (Figure 4A, Supplementary Table S6). Initially, the univariate Cox regression analysis revealed nine genes associated with the HCC prognosis (Supplementary Table S7). Finally, LASSO regression analysis identified eight prognostic signature genes: PDZ and LIM domain 3 (PDLIM3), peptidylglycine alpha-amidating monooxygenase (PAM), PDZ and LIM domain 7 (PDLIM7), fascin actin-bundling protein 1 (FSCN1), dihydropyrimidinase like 2 (DPYSL2), AT-rich interaction domain 5B (ARID5B), galectin 3 (LGALS3), and Kruppel-like factor 2 (KLF2) (Figures 4B, C). The prognostic genes were enriched significantly in 403 terms according to the GO database, the top-10 of which were displayed in a bubble plot (Figure 4D). The biological process (BP) category included “multicellular organism growth” and “T cell activation via T cell receptor contact with antigen bound to MHC molecule on antigen presenting cell”. The cell component (CC) category included “stress fiber” and “contractile actin filament bundle”. The molecular function (MF) category included “muscle alpha-actinin binding” and “alpha-actinin binding”. Moreover, TCGA-LIHC samples were consistently clustered into different clusters according to the expression of prognostic genes. It can be seen that the area under the cumulative density function (CDF) curve increased significantly when k ≤ 4 (Figure 4E), but the area under the CDF curve did not increase significantly when k≥5. And when k=5, the effect of consensus clustering was not good (Supplementary Figure S2). Therefore, HCC patients were divided into 4 clusters (Figure 4F). Differences in gene expression and clinicopathological characteristics among these four clusters were demonstrated with a heatmap (Figure 4G). Most importantly, there was a significant difference in survival between the four clusters (P=0.002) (Figure 4H), which initially demonstrated the prognostic value of these eight genes.

According to the coefficients (Table 1) and expression of prognostic genes, the risk score of each sample in TCGA-HCC was calculated. Then, HCC patients in the training set were divided into low- and high-risk groups based on the optimal cutoff of risk score (0.126) (Figure 5A). Overall survival (OS) (Figure 5B), DFS, PFS, and DSS (Supplementary Figures S3A–C) were longer in the low-risk group than in the high-risk group, which indicated that patients in the low-risk group had a better overall prognosis. To evaluate the performance of the risk model, ROC curves were plotted, and the area under the ROC curve (AUC) at 1, 3, and 5 years was 0.728, 0.689, and 0.663, respectively (Figure 5C). The results for univariate and multivariate analyses (Figures 5D, E) and the Concordance index (C-index) (Figure 5F) showed that the risk score was: (i) an independent factor affecting survival; (ii) a superior prognostic predictor than other indicators. Moreover, GSEA showed that “epithelial cell signaling in Helicobacter pylori infection”, “regulation of actin cytoskeleton”, “p53 signaling pathway”, “focal adhesion” and “MAPK signaling pathway” were enriched significantly in the high-risk group, whereas “fatty acid metabolism”, “glycine, serine and threonine metabolism”, “primary bile acid biosynthesis”, “tryptophan metabolism” and “valine, leucine and isoleucine degradation” were enriched significantly in the low-risk group (Figure 5G).

To verify the reliability of the prognostic signature, we further validated it in the test set. Samples were grouped in the same way as in the training set (Figure 6A). Patients in the high-risk group had a worse prognosis than those in the low-risk group (Figure 6B), and had an AUC of 0.701, 0.677, and 0.653 at 1, 3, and 5 years in the test set, respectively (Figure 6C). These results validated the reliability of the M2-like TAM-related prognostic signature in predicting the prognosis of HCC patients. Based on the risk score of our prognostic signature and other clinicopathological indicators of patients, we constructed a nomogram to make a more comprehensive prediction of patient survival (Figure 6D). Moreover, the results of the calibration curve and ROC curve of the nomogram showed a reliable performance, with an AUC of 0.764, 0.730, and 0.737 at 1, 3, and 5 years, respectively (Figures 6E, F).

In addition to significant differences in survival between high- and low-risk groups, they also differed in their clinicopathological characteristics. Figure 7A shows a heatmap of the clinicopathological characteristics and expression of signature-related genes in high- and low-risk groups. There were no significant differences between high- and low-risk groups in terms of sex and age distribution, whereas there were significant differences in the depth of tumor infiltration (T stage) and tumor grade, with a significantly higher proportion of patients of grade 3 and T3–4 in the high-risk group than in the low-risk group (Figure 7B). Survival analyses of HCC patients divided into different subgroups according to their clinicopathological indicators showed that the survival outcome of patients in the high-risk group was worse than that of the low-risk group, whether grouped by sex, age, grade, or T stage (Figures 7C–J).

Based on risk grouping, in terms of immune cells, we discovered that the content of macrophages and T_regs was higher in the high-risk group than that in the low-risk group, whereas the content of B cells, mast cells, NK cells, plasmacytoid dendritic cells (pDCs), and helper T cells was lower than that in the low-risk group (Figure 8A). In terms of immune functions, the high-risk group was lower than the low-risk group in terms of cytolytic activity and a type-II interferon response, but higher than the low-risk group for major histocompatibility class (MHC) class-I (Figure 8B). With regard to immunotherapy, TIDE scores were higher in the low-risk group than those in the high-risk group (Figure 8C), suggesting that patients in the low-risk group were more likely to experience immune evasion, and that immunotherapy may be less efficacious. There was no significant difference in the scoring of several immunotherapy treatments between patients in high- and low-risk groups (Supplementary Figures 1A–D). In terms of expression of the genes associated with immune checkpoints, many differentially expressed genes between the two groups were documented (Figure 8D), such as CD44, CD86, and CD276, which showed significantly higher expression in the high-risk group. This finding offers the possibility of discovering new targets for immunotherapy.

To further investigate the clinical use of prognostic genes, we employed the CellMiner database to explore the relationship between prognostic genes and drug sensitivity. PAM was correlated significantly and positively with simvastatin sensitivity (correlation = 0.442, P<0.001) and LGALS3 was correlated significantly and positively with ARRY-162 sensitivity (correlation = 0.414, P<0.001) (Figure 9A). Patients in the high-risk group had a significantly worse prognosis, so we predicted 10 drugs with higher sensitivity in the high-risk group: epothilone B, A-443654, BEZ235, BI-2536, BMS-75480, CGP-6047, foretinib, GSK212645, JW-7-52-1, and VX-680. Epothilone B had the lowest IC₅₀ (Figures 9B–K).

After obtaining the M2-like TAM-related biomarkers and constructing related prognostic signature, we further analyzed the expression of signature genes in TCGA-LIHC samples. Figure 10A showed that the RNA expression levels of PDLIM3, PAM, PDLIM7, FSCN1, and LGALS3 in tumor samples were significantly upregulated. Moreover, in the samples we obtained from HCC patients, the RNA expression levels of these 5 genes (Figures 10B–F) were also significantly higher in tumor tissues than in adjacent normal tissues, perhaps suggesting that these genes play a role in the progression of HCC.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/. 12