Identification of a Macrophage Marker Gene Signature to Evaluate Immune Infiltration and Therapeutic Response in Hepatocellular Carcinoma Based on Single-cell and Bulk RNA Sequencing Data

doi:10.21203/rs.3.rs-3561035/v1

Download PDF

Research Article

Identification of a Macrophage Marker Gene Signature to Evaluate Immune Infiltration and Therapeutic Response in Hepatocellular Carcinoma Based on Single-cell and Bulk RNA Sequencing Data

https://doi.org/10.21203/rs.3.rs-3561035/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Only a minority of hepatocellular carcinoma (HCC) patients can benefit from systemic regimens. Macrophages, which abundantly infiltrate in HCC, could mediate tumour microenvironment remodelling and immune escape, proving to be powerful weapons in combating HCC. Thus, a deeper understanding of macrophages is necessary for improving existing antitumour treatments.

Methods

By leveraging a series of bioinformatic and machine learning approaches, we comprehensively explored the transcriptomic landscape of human HCCs and the role of macrophages related genes in HCC from multiple single-cell and bulk RNA sequencing datasets. The uniform manifold approximation and projection (UMAP) was utilized for dimensionality reduction and the FindAllMarkers function was used to detecting for differentially expressed genes (DEGs). Unsupervised clustering was performed using the “ConsensusClusterPlus” to cluster the macrophage marker genes (MMGs). Further, Gene set variation analysis (GSVA) and Function Enrichment Analysis to elucidate the functional difference of the MMGs associated clusters. Subsequently, identified 12 hub genes, based on which component analysis algorithm (PCA) was employed to construct Macrosig score. Meanwhile, the prognosis, biological characteristics, mutation profile, TME cell-infiltrating and drug response of different Macrosig score were analysised.

Results

Our studies identified 13 macrophage marker genes (MMGs) based on 574 HCC samples, based on which three MMG-associated clusters(A,B,C) were defined. Patients in cluster B were significantly associated with favorable overall survival, clinicopathological features and a lower immune infiltration score. On this basis, 12 hub genes were identified among these clusters, subsequently a scoring system by principal component analysis was constructed terming the Macrosig score. Importantly, patients with low Macrosig score, characterized by higher immune infiltration, increased mutation frequency and higher immune checkpoint expression including CTLA-4, LAG3, PDCD1 and TIGIT, exhibiting enhanced efficacy of immunotherapy as expected when validated in external database. Moreover, for molecular targeted therapy, low Macrosig score indicates increased sensitivity to AZD.2281, A.443654, ABT.263, ABT.888, AG.014699 and ATRA, while high Macrosig score were more likely to be sensitive to AZD6482, AKT inhibitor VIII, AS601245, AZ628, AZD.0530 and AZD6244.

Conclusions

In summary, a novel scoring system based on macrophage subclusters was constructed, thereby guiding more effective prognostic evaluation and tailored therapeutic regimens of HCC patients.

Hepatocellular carcinoma(HCC)

Single-cell RNA sequencing

Macrophages

Scoring system

Drug therapy

Hepatocellular carcinoma (HCC) is the most common primary malignancy of the liver and the fourth leading cause of cancer-related mortality worldwide[1]. Over the past two decades, various studies have examined the clinical management of HCC, and novel systemic treatment options have emerged. In particular, the PD-L1 inhibitor atezolizumab plus the antiangiogenic agent bevacizumab as well as the CTLA-4 inhibitor tremelimumab-actl plus the PD-L1 inhibitor durvalumab have been approved as preferred first-line systemic therapy regimens for HCC, indicating that there is increasing support the use of immune checkpoint inhibitor (ICI)-related regimens[2–4]. However, despite the considerable advances that have occurred in recent years, several questions remain unanswered: although monotherapy with PD-1- or PD-L1-related agents is generally well tolerated, the risk of immune-related adverse events is increased with combination regimens. Furthermore, the number of HCC patients who have minimal or no clinical benefit from such agents is still large. The greatest focus has been the expression of immune checkpoints on tumour cells, and although positive expression of such checkpoints can be used to identify populations who are likely to benefit clinically, testing alone is not sufficient for HCC patient selection[5]. Therefore, reliable predictors need to be identified to optimize patient benefit, minimize the risk of toxicities and guide combination approaches in HCC.

As a relatively immunotherapy-resistant malignancy characterized by high cellular heterogeneity, HCC tumours are composed of a complex tumour microenvironment (TME) consisting of cellular (tumour-infiltrating immune cells and stromal cells), chemical (cytokines and chemokines), and extracellular matrix components[6, 7]. The TME is thought to play a crucial role in HCC, and the multifaceted interactions between tumour cells and surrounding immune components in the TME can affect tumour progression, inhibit immune surveillance, and induce immune escape and tolerance by inducing the secretion of various molecules[8]. In particular, substantial clinical and experimental evidence indicates that there is a substantial level of macrophage infiltration in the TME in HCC tumours. Moreover, as fundamental immune cells in innate and adaptive immunity, macrophages have a dual role in liver cancer, and they are easily driven towards immunosuppressive tumour-associated macrophages (TAMs) upon interactions with TME components, which are strongly associated with patient survival and correlate with drug resistance[9]. In mouse models of cancer, macrophages promote cancer initiation and progression by stimulating angiogenesis, increasing the migration, invasion and intravasation of tumour cells and suppressing antitumour immunity[10]. Other studies have also suggested that macrophages are involved in tumour antigen recognition disorders, the recruitment and function of immunosuppressive cells, the secretion of immunosuppressive cytokines, immune checkpoint crosstalk, and the formation of immune privileged sites and can directly suppress T-cell responses via PD-L1 in HCC[11, 12]. Notably, although macrophages express canonical markers, studies have indicated considerable transcriptomic diversity between macrophage populations even within a tissue. Given the dominant roles of macrophages in immunity, a deeper understanding of the immune cell infiltration characteristics mediated by macrophages will improve tumour immunotherapy.

Furthermore, advances in next-generation sequencing have led to an in-depth understanding of the genomics and epigenomics of HCC. Likewise, recent advances in single-cell RNA sequencing (scRNA-seq) have overcome some of the technical hurdles in the investigation of cellular heterogeneity among complex tissues such as carcinoma tissues, contributing to the identification of different subclasses of immune environments that have an influence on tumour initiation and the response and therapy. In addition, the availability of genomic datasets from public databases has made it much easier to study numerous transcriptomics and genomic profiles, providing a means to comprehensively analyse the interactions between macrophages and immune regulatory factors to predict immunotherapy response and guide treatment[13, 14].

Herein, we take advantage of the unique multidimensional capacity of scRNA-seq to explore the landscapes of immune cell infiltration in human HCCs, comprehensively delineate the multifaceted nature of macrophages, and evaluate the association between macrophage subclusters and infiltrating immune cells in the TME by integrating TCGA-liver hepatocellular carcinoma (LIHC), GEO, HCCDB and TISCH transcriptomic and genomic datasets. In this way, we first identified 13 macrophage marker genes (MMGs), based on which three distinct MMG-associated clusters were classified. Moreover, we further constructed the Macrosig scoring system to predict prognosis and response to immunotherapy and targeted drugs, which can be used as a tool to improve HCC therapy.

Raw dataset acquisition

A total of five scRNA-seq datasets were acquired from TISCH (http://tisch.comp-genomics.org/) in this study. The bulk RNA-seq profile of the LIHC cohort was employed as a training cohort for model development. Publicly available gene expression data and complete clinical annotations from patients with LIHC were retrospectively obtained from the HCCDB (http://lifeome.net/database/hccdb/home.html), and the TCGA-LIHC data were obtained from the UCSC-XENA website (https://xenabrowser.net/datapages/). Patients without overall survival (OS) information, samples derived from borderline tumours, benign tissues or normal liver tissues and GEO samples from TCGA were excluded from this work, and a total of 574 samples with 17661 genes were included in the subsequent analyses. Furthermore, two large cohorts from the GEO database were selected as external validation cohorts. The GSE176307 dataset, which included the detailed clinical information and expression matrix of 103 patients with metastatic urothelial cancer (mUC) receiving immune checkpoint blockade (ICB) agents, was also obtained from the GEO database for further analysis[15]. In addition, data from the phase 3 JAVELIN Renal 101 trial (NCT02684006) included the complete clinical and expression profiles of advanced renal cell carcinoma (aRCC) patients receiving avelumab (anti-PD-L1)[16] and were acquired via R software.

Simple nucleotide variation (SNV) data from TCGA were also downloaded for tumour mutation burden (TMB) calculation. All RNA-seq data were converted to transcripts per million (TPM) format and further log2 transformed. The Limma and SVA algorithms were applied for raw data batch normalization[17, 18].

Bulk RNA-seq and scRNAseq data processing and analysis

According to the standardized pipeline, the 10X scRNA-seq data were processed through R software (version 4.2.2) via the “Seurat” package[19]. Quality control (QC) was performed on the raw matrix to obtain a high-quality scRNA-seq expression matrix for further analysis by filtering out low-quality cells according to the following criteria: 1) genes that were expressed in at least five single cells and cells that expressed more than 300 genes were selected to create a Seurat object; 2) only cells that expressed more than 300 genes and less than 7500 genes were included; 3) the percentage of mitochondria of each cell was calculated, and cells with more than 4% mitochondrial genes were regarded as low-quality cells and were excluded from downstream analysis; and 4) cells with a percentage of ribosomal genes more than 3% and percentage of haemoglobin less than 0.1% were included. After filtering, 22,143 genes and 59,509 cells were included. In addition, the “harmony” method in the “ScaleData” function was used to normalize the scRNA-seq data. Subsequently, the “RunPCA” function in the “Seurat” package was utilized for principal component analysis (PCA). After that, uniform manifold approximation and projection (UMAP) was used for dimensionality reduction and cluster identification. Furthermore, the “FindAllMarkers” function in the “Seurat” package was exploited for cell cluster identification, with the parameter “resolution” set as 0.2. Then, the significant differentially expressed genes (DEGs) of each cluster were identified by calculating the log2 [fold change (FC)] value and the adjusted p value. DEGs with |log₂FC|≥1 and adjusted p value < 0.05 were considered MMGs of each cluster. Finally, the “SingleR” package was employed for automatic cluster annotation to identify the cell types based on information from the Human Primary Cell Atlas[20].

Identification of MMG-associated clusters and estimation of immune cell infiltration

Unsupervised clustering was performed using the “ConsensusClusterPlus” package in R to cluster the MMGs[21] with 1,000 repetitions to ensure the stability of classification. Three MMG-associated clusters were identified. The Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data (ESTIMATE) algorithm, which can be used to estimate tumour cellularity as well as the tumour purity according to unique transcriptional profile characteristics, was used to calculate the ImmuneScore, StromalScore and ESTIMATE score and predict the levels of infiltrating immune cells and tumour purity[22]; furthermore, their correlations with specific cell types were analysed by correlation analysis, with p value < 0.05 being statistically significant. Moreover, the Wilcox test was performed to investigate the relationship between MMG subclusters and common HCC clinicopathological characteristics, including age, sex, clinical stage, T stage and survival outcome in the TCGA cohort.

GSVA and functional enrichment analysis

To elucidate the functional differences in the three MMG-associated clusters obtained from the former analysis, we performed Hallmarks of Cancer[23], Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome enrichment analyses via data from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/) by using the “clusterProfiler” R package to analyse biological functions and comprehensively scored each gene set by GSVA to evaluate differences in biological functions between different samples at the gene and pathway levels[24, 25]. We set adjusted p value < 0.05 and |logFC|>1 as cut-off criteria, and 761 DEGs were identified. A heatmap was used to illustrate differences in pathway activity between two related subclusters.

Analyses of the hub genes and mutation landscapes of the MMG clusters

To clarify the prognostic value of the identified MMG clusters in HCC, we performed univariable Cox regression analysis using the 761 DEGs, and 12 hub genes of the 761 DEGs with p < 0.000001 were identified among the three MMG clusters: G6PD, MYBL2, LPCAT1, PKIB, ATP1B3, PON1, AFM, SLC1A5, SPP1, SLC22A1, ADH1C, and SMOX; these hub genes also demonstrated promising prognostic power in the Cox regression analysis. Gene mutation data of HCC patients were downloaded from TCGA, and TMB was calculated using the “maftools” package. We further comprehensively analysed the 12 hub genes using GSCA (http://bioinfo.life.hust.edu.cn/GSCA/#/), including analysing somatic mutation and copy number variation (CNV) profiles downloaded from the Genomic Data Commons. Somatic mutation data calculated by the R package “maftools” were sorted in mutation annotation format (MAF) format.

Construction and validation of a the prognosis-related Macrosig scoring system

To quantify the MMG-related profiles of individual tumours, a scoring system was developed, termed the Macrosig scoring system. PCA was employed to establish the Macrosig scoring system based on the expression of the above 12 MMG-associated prognostic hub genes, and principal components PC1 and PC2 were identified. The sum of PC1 and PC2 served as the final Macrosig score. The testing and validation samples were divided into low- and high-score groups, each of which was subjected to immune infiltration analysis, and Kaplan‒Meier survival curves were used to visualize the survival difference between the groups with high and low Macrosig scores using the “survminer” package in R software.

Immunotherapy and targeted drug response prediction

To explore the correlation of the Macrosig score with immunotherapy response, a boxplot was used to display differences in the expression of four immune checkpoints between the two Macrosig score groups. Then, two independent datasets, GSE176307 and NCT02684006, with complete clinical information from patients who received immunotherapy were downloaded and analysed to determine the utility of the Macrosig score for predicting immunotherapy response. Next, to identify differences in the effects of targeted drugs between patients in the two groups, the half maximal inhibitory concentration (IC50) values of specific chemotherapeutic drugs commonly used for HCC treatment were estimated using the package “pRRophetic” based on data from Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/)

Statistical analysis

All statistical analyses and data visualization were conducted using R software version 4.1.1 for Windows 64.0 (http://www.R-project.org) in this study. The nonparametric Wilcoxon t test was used to investigate differences in continuous categorical variables between the two groups. Comparisons of proportions were performed with the chi-square test or Fisher’s exact test. Cox regression analysis was used to identify prognostic genes. Survival analysis was performed by the Kaplan‒Meier method. Unless otherwise noted, p < 0.05 was set as the significance threshold. The Benjamini–Hochberg method was used to adjust the p value for multiple tests using the R function “p.adjust”.

The landscape of immune cell infiltration in HCC as assessed by scRNA-seq analysis

Based on scRNA-seq data from GSE140228[26], we obtained gene expression profiles of 59509 high-quality cells from sixteen patients who were pathologically diagnosed with liver cancer, including thirteen males and three females, and these cells were identified for subsequent analysis after strict quality control. Fourteen cell clusters were identified by dimensionality reduction analysis with UMAP. Subsequently, the cell identity of each cluster was annotated using a reference dataset from the Human Primary Cell Atlas, and six cell types were identified: natural killer (NK) cells, T cells, mononuclear macrophages, B cells, plasma cells and endothelial cells (Fig. 1A). Next, the proportions and relationships of different cell types in each cluster were identified, as displayed in Fig. 1B; T cells, NK cells and mononuclear macrophages were the three most predominant cell types in HCC. The FindAllMarkers function was used to detect DEGs in the different cell types. The top 5 genes with high or low expression ranked according to the log₂FC value for each of the six types of cells are shown in Fig. 1C, and all the differentially expressed marker genes are listed in Additional file 5. Table. S1.

We further conducted GSVA using the hallmark gene sets as reference gene sets derived from MSigDB to determine the enrichment scores of the different immune cell clusters to provide insights into the different biological functions of the groups. Excitingly, we found that mononuclear macrophages showed enrichment of the following hallmarks: inflammatory response, complement, TNF-α signalling via NF-κB, oxidative phosphorylation, adipogenesis and mTORC1 signalling; thus, these cells had a promising antitumour immunophenotype (Fig. 2). Therefore, in this study, we focused on the identification of a mononuclear macrophage marker genes (MMG) signature and further investigated the relationships between different immune subtypes and molecular subtypes.

Identification of MMG expression profiles in HCC

To unravel the underlying biological characteristics of macrophage immunophenotypes, we identified a total of 17661 genes of 574 samples from the TCGA expression profile dataset, with 371 tumour samples, and the HCCDB18 dataset, with 203 tumour samples. The limma and SVA algorithms were applied to the data for batch correction, and two datasets were separated in the PCA plots because of heterogeneity, which was eliminated after normalization (Fig. 3A). Next, we selected the top 50 MMGs with high expression in macrophages identified by the abovementioned scRNA-seq analysis, among which 48 genes were expressed in the combined data. 32 of the 48 MMGs showed a significant effect on patient survival outcomes in the combined dataset (Additional file 1. Fig. S1). Ultimately, we obtained 13 MMGs in HCC according to cut-off criteria of |logFC|>1 and adjusted p value<0.05 in the univariate Cox regression analysis (Additional file 6. Table. S2) : S100A9, CD14, GRN, NPC2, CD68, IER3, MAFB, HLA-DRB5, CPVL, IFI30, PLAUR, S100A8, and FCER1G; furthermore, the comprehensive landscape of MMG interactions, regulatory connections, and prognostic value in patients with HCC is shown in Fig. 3B. We next analysed the expression distribution of these 13 MMGs, and as expected, they were mainly specifically expressed in the macrophage cell cluster (Fig. 3C). In addition, we combined multiple single-cell liver cancer datasets from TISCH to verify the above conclusion and confirmed that the 13 identified MMGs were generally highly expressed in macrophages (Additional file 2. Fig. S2).

Three MMG-associated clusters and their functional differences

To further explore the expression profiles of tumour-infiltrating MMGs in HCC, we applied an unsupervised consensus clustering algorithm to categorize patients with HCC based on the expression of the 13 MMGs. Our results showed that k=3 was the optimal selection value for sorting the entire cohort into 3 clusters: A (n=145), B (n=194) and C (n=227) (Fig. 4A). Furthermore, we explored the prognostic implications of macrophage gene clusters by integrating them with survival information. The analysis was performed by using Kaplan‒Meier plotter, and we found that patients in gene cluster B had a better prognosis, whereas patients in gene clusters A and C had unfavourable outcomes (log rank test, p=0.028; Fig. 4B). In addition, we found that the expression levels of MMGs were mostly lower in cluster B than in clusters A and C, except for that of CD14 (Fig. 4C). Considering the role of macrophage subtypes in HCC, we also compared the clinical features of the 3 clusters of macrophages. The relationships between HCC patient clinical factors and macrophage subtypes were delineated in a heatmap (Fig. 4D).

To identify differentially enriched pathways among the 3 macrophage clusters in the HCC cohort, we downloaded the Hallmark, KEGG and Reactome pathway sets from MSigDB and carried out GSVA to identify the enriched pathways of these 3 clusters. The enrichment heatmap illustrated that cluster B had significant enrichment of inflammatory and antitumour hallmarks (including PI3K-AKT-mTOR signalling, IFN-γ response, IFN-α response, IL-6-JAK-STAT3 signalling, IL-2-STAT5 signalling, Kras signalling, TNF-α signalling via NF-κB and inflammatory response). Cluster B had significant enrichment of KEGG pathways related to immune infiltration (including NOD-like receptor signalling pathway, Toll-like receptor signalling pathway, T-cell receptor signalling pathway, B-cell receptor signalling pathway, cytokine‒cytokine receptor interaction and chemokine signalling pathway). For Reactome pathways, cluster B had significant enrichment of inflammation and immune activity (including PD-1 signalling, DAP12 interactions, IL-10 signalling, and signalling by interleukins and inflammasomes) (Fig. 5).

Estimation of TME infiltration and functional differences in distinct macrophage clusters

Macrophages play a vital role in antitumour immunity. Though a close relationship between macrophage subclusters and immune activity was identified, macrophages and the immune response are processes involving multiple steps and molecules that act in a highly united and cooperative manner. We further explored the relationship of the macrophage subclusters with immune cell infiltration in HCC patients. PCA was used for dimensionality reduction (Fig. 6A), and by using the ESTIMATE algorithm, we found that cluster B had lower StromalScores, ImmuneScores and ESTIMATE scores than clusters A and C (p<0.001) (Fig. 6B). To investigate the role of MMGs in the TME of HCC, we assessed the correlations between the 3 macrophage clusters and 23 human immune cell subsets using single-sample gene set enrichment analysis (ssGSEA). Notably, we observed significant differences in the infiltration of most immune cells between the three clusters; for example, cluster A was found to have more immune cell infiltration, but cluster B had less immune infiltration (Fig. 6C). Next, we performed differential gene analysis for the three clusters, and DEGs among the three clusters were visualized in volcano plots. As shown in Fig. 6D, the number of DEGs between cluster B and cluster A was the largest, which was consistent with the above phenotypic results. On this basis, the gene sets were filtered using the criteria |logFC|> 1 and p<0.05, and 761 DEGs were acquired. The clusterProfiler R package was adopted for GO analysis. The 761 DEGs were mainly enriched in the biological process (BP) terms cell‒cell adhesion, leukocyte activation, and cell chemotaxis; the cellular component (CC) terms collagen-containing extracellular matrix, external side of plasma membrane and vesicle lumen; and the for molecular function (MF) terms glycosaminoglycan binding, extracellular matrix structural constituent and immune receptor activity (Fig. 6E). Further KEGG pathway analysis found that the DEGs were mainly enriched in complement and coagulation cascades, phagosomes, drug metabolism (cytochrome P450), and metabolism of xenobiotics by cytochrome p450, which may imply that the macrophage subclusters have effects on antitumour drug metabolism (Fig. 6F).

Hub gene subtypes based on macrophage clusters and construction of the Macrosig scoring system

Multiple studies have shown that the immune system is related to therapy outcomes, and protumour or antitumour activity can exist [27, 28]. To clarify the prognostic value of the identified macrophage clusters in HCC, we performed a univariable Cox regression analysis using the 761 DEGs, and 12 hub genes of the 761 DEGs with p<0.000001 were identified among the three macrophage clusters (G6PD, MYBL2, LPCAT1, PKIB, ATP1B3, PON1, AFM, SLC1A5, SPP1, SLC22A1, ADH1C, and SMOX), which also demonstrated noticeable prognostic power in the Cox regression analysis (Fig. 7A, Additional file 7. Table S3). In addition, we also conducted a comprehensive analysis of 12 genes using the GSCA datasets, and they were all differentially expressed in the HCC TCGA dataset; we also analysed deleterious mutations, mutation distribution, CNV percentage, and methylation differences (Additional file 3. Fig. S3).

Given the crucial regulatory function of macrophages in the evaluation of prognosis and TME immune landscapes, we further constructed a scoring system with PCA (PC1+PC2) based on the 12 hub DEGs. We termed it the Macrosig scoring system. We divided patients into low score and high score groups to quantify macrophage regulation and the immune microenvironment in individual HCC samples. We found that a higher Macrosig score indicated a longer survival time (p<0.001) (Fig. 7B). Moreover, according to the tumour-infiltrating immunocyte analysis, we determined the levels of 23 types of infiltrating immune cells in patients with different Macrosig scores; we found that the Macrosig score was negatively correlated with the majority of infiltrating immune cells (Fig. 7C). These results preliminarily suggest that the higher the Macrosig score is, the less immune cell infiltration and the better the prognosis of HCC patients.

Correlation between the Macrosig score, immunological characteristics and clinical features

In terms of clinicopathological features, the Macrosig score was also correlated with patient stage and survival status. As shown in Fig. 7D, significant differences in the Macrosig score were observed among different HCC stages; patients with a low Macrosig score (59% stage I+II and 41% stage III+IV) presented more severe stages of disease than those with a high Macrosig score (77% stage I+II and 23% stage III+IV) (p<0.001). More importantly, we also observed that the high Macrosig score group (80% alive and 20% dead) presented a prominent survival advantage over the low Macrosig score group (55% alive and 45% dead) (p<0.001, Fig. 7E), indicating that the Macrosig score has potential prognostic value for HCC patients. Furthermore, we investigated the correlation between the Macrosig score and immune checkpoints. In line with the above negative correlation of the Macrosig score with the levels of most infiltrating immune cells, the results showed that patients with a low Macrosig score had significantly higher expression of most immune checkpoints, including CTLA-4, LAG3, PDCD1 and TIGIT (Fig. 7F), than patients with a high score, which suggested that HCC patients in the low Macrosig score group were more likely to be responsive to immunotherapy.

Increasing evidence has demonstrated that somatic mutation patterns are associated with responsiveness to immunotherapy. We thus assessed the distribution of somatic variants in HCC driver genes between the low and high Macrosig score subgroups. The HCC driver genes were identified by using the maftools R package[29], and the top 10 driver genes with the highest alteration frequency in the high (altered in 157 of 228 samples) and low (altered in 95 of 124 samples) Macrosig score groups were further analysed. The mutation annotation files revealed that the alteration frequencies of TP53, NFASC, EPHA4, DMD, CPAMD8, NLRP2, ADAMTSL3, ADARB2, IQGAP3, LPIN2 and SLC6A11 were significantly higher in the low score group than in the high score group, but this pattern did not hold for DYNC2H1 (Fig. 8A-C). These results might provide novel ideas for investigating the mechanism of gene mutation in immune checkpoint blockade therapy.

The relationship of the Macrosig score with cytokines

Chemokines, interleukins and interferons are all cytokines; they play an important role in the composition of the TME, coordinate immune cell trafficking, exert multifaceted roles in the tumour immune response and tumour biology, and have a strong influence on patient prognosis and response to therapy; therefore, the cytokine network has emerged as a potential immunotherapy target[30-32]. TAMs generally include two subtypes: M1 macrophages play an antitumour immunity role, while M2 macrophages promote tumour growth, invasion and migration and inhibit tumour immunity by secreting anti-inflammatory cytokines. However, whether the Macrosig score has an effect on cytokine expression is unknown. As such, we next examined the expression of various chemokines in the high and low Macrosig score groups and observed that the expression of most immune activity-related cytokines was significantly different between the groups with different Macrosig score, with the high score group expressing decreased levels of cytokines (Fig. 9A). In the subsequent analysis, the correlations between the Macrosig score and 50 hallmark terms were calculated, and the results were visualized via ggplot2 (Fig. 9B). The Macrosig score was positively correlated with xenobiotic metabolism, bile acid metabolism, fatty acid metabolism, and so on, and negatively correlated with MYC targets V1, G2M checkpoint, E2F targets, unfolded protein response, unfolded protein response, DNA repair, and so on.

Macrosig score-related treatment strategy for HCC

Immune checkpoint blockade therapy, which is used to block T-cell inhibitory molecules in cancer treatment, has shown promising outcomes and the potential to improve disease conditions in advanced cancers, but it is not effective in all patients[33, 34]. The higher the Macrosig score was, the lower the immune infiltration, immune checkpoint expression, cytokine expression and frequency of gene mutation. Therefore, we hypothesized that the high score group would be insensitive to immunotherapy. To test this, we used the immunotherapy database to verify the efficacy of the Macrosig score in predicting immunotherapy response. By employing a real-world dataset derived from patients with metastatic urothelial cancer treated with immune checkpoint inhibition (GSE176307)[15], we discovered that the survival probability was significantly higher in the low Macrosig score group (p=0.046); moreover, the complete response (CR) or partial response (PR) rate was also higher than that in the high Macrosig score group (29% CR+PR vs. 12% CR+PR) (Fig. 10A, B). In addition, based on data from the phase 3 JAVELIN Renal 101 trial (n = 886; NCT02684006) that aimed to evaluate the efficacy of avelumab (anti-PDL1 antibody)[16], we further revealed that patients with low Macrosig scores exhibited longer survival times and lower disease progression rates than those with high Macrosig scores (p=0.0049) (Fig. 10C, D); all of these results demonstrated that the Macrosig score can play a crucial role in predicting the immune response.

With regard to targeted immunotherapy, we comprehensively evaluated the sensitivities of targeted drugs in different Macrosig score groups using pRRophetic, with a higher IC50 indicating less treatment sensitivity. The results showed that patients in the low score group had lower IC50 values for chemotherapeutics, including AZD.2281 (olaparib), A.443654 (an Akt inhibitor), ABT.263 (navitoclax), ABT.888 (a PARP1 inhibitor), AG.014699 (rucaparib) and ATRA; that is, patients in the low Macrosig score group were more likely to benefit from these drugs. The IC50 values of AZD6482 (a PI3Kβ inhibitor), AKT inhibitor VIII, AS601245 (a JNK inhibitor), AZ628 (a Raf kinase inhibitor), AZD.0530 (saracatinib), and AZD6244 (selumetinib) were significantly lower in the patients in the high Macrosig score group than in those in the low Macrosig score group (Fig.6E), suggesting that high Macrosig score patients are more likely to be sensitive to these treatments. These results indicate that the Macrosig score can be used to predict sensitivity to chemotherapeutic drugs and predict patient treatment benefit. In addition, we used the RegNetwork database to predict the upstream miRNAs and transcription factors of the above 12 genes (Additional file 4. Fig. S4) to lay a foundation for further mechanistic studies.

HCC is among the tumour types that are resistant to immune checkpoint inhibition, with only a minority of HCC patients benefiting. Finding reliable predictors to optimize patient benefit, minimize the risk of toxicities and guide combination approaches for HCC patients remains a challenge. Mounting evidence has demonstrated that macrophages are highly enriched in the tumour immune microenvironment and play an indispensable role in inducing immunosuppression in liver malignancies[35, 36]. Previous studies have verified that macrophages act as key orchestrators associated with cancer progression and distant metastasis, and their abundance is correlated with clinical prognosis; for example, macrophages can promote the expression of PD-L1 on HCC cells through cytokines or directly suppress T-cell responses via PD-L1, and cytokines secreted by macrophages can affect TAM-mediated phagocytosis via the STAT3 signalling pathway, resulting in changes in overall survival and recurrence-free survival[11, 12, 37]. However, there is still significant room for improvement in macrophage-mediated immune modulation, and a comprehensive molecular analysis of macrophages in HCC is relatively lacking[38, 39].

On the basis of existing research, herein, we carried out a comprehensive bioinformatic analysis of the transcriptomic landscape of human HCCs based on scRNA-seq and bulk RNA sequencing data. In the present study, we obtained 13 MMGs (S100A9, CD14, GRN, NPC2, CD68, IER3, MAFB, HLA-DRB5, CPVL, IFI30, PLAUR, S100A8, FCER1G), all of which are mainly specifically expressed in the macrophage cell cluster and correlated with the prognosis of HCC patients; this finding may provide insights into the molecular mechanisms for further clinical studies to improve antitumour therapy efficacy in HCC. Next, we identified 3 distinct MMG-associated clusters (clusters A, B and C) and evaluated the association between these subclusters and infiltrating immune cells in the TME. We observed that cluster B patients had more favourable outcomes than those in clusters A and C. As expected, cluster B was significantly characterized by enrichment of inflammatory and immune-related biological behaviours and pathway terms, including the NOD-like receptor signalling pathway, Toll-like receptor signalling pathway, T-cell receptor signalling pathway, B-cell receptor signalling pathway, cytokine‒cytokine receptor interaction and chemokine signalling pathway, PD-1 signalling, DAP12 interactions, and signalling by interleukins and inflammasomes. We also identified that cluster B had a lower StromalScore, ImmuneScore and ESTIMATE score, in line with a previous study showing that a high immune infiltration score was associated with a poor prognosis, confirming the reliability of our classification of immune phenotypes for the different clusters[40].

Moreover, 761 DEGs were identified among the three clusters, and functional enrichment analysis showed that DEGs were mainly enriched in the terms complement and coagulation cascades, phagosome, drug metabolism (cytochrome P450), and metabolism of xenobiotics by cytochrome p450, implying that the macrophage subclusters may affect antitumour drug metabolism. The 12 hub genes (G6PD, MYBL2, LPCAT1, PKIB, ATP1B3, PON1, AFM, SLC1A5, SPP1, SLC22A1, ADH1C, SMOX) were further assessed. Given the crucial regulatory function of macrophages in prognostic evaluation and the potential effects of macrophage subclusters on antitumour drug metabolism identified in this study, a comprehensive evaluation of macrophage subtypes will enhance our understanding of TME cell infiltration. Moreover, we established a reliable and effective system termed as the Macrosig scoring system to evaluate macrophages in individual HCC patients. Notably, we found this system to be an efficient prognostic tool that can guide precision immunotherapy and targeted drug therapy. Our research lays a foundation for a deeper understanding of patient antitumour immune response and provides guidance for more individualized and effective combination therapy regimens related to the Macrosig score.

However, there are still some inevitable limitations and shortcomings of our research: scRNA-seq and bulk RNA-seq data were combined to explore the landscapes of immune cell infiltration and comprehensively delineate the multifaceted nature of macrophages in human HCC. In spite of the accuracy of transcriptome data is relatively high, and the sample size of the scRNA-seq data was limited. The TME is shaped by various cell types and factors, and there is a complex network among different components. Therefore, the immune cell infiltration and immune function status estimated by algorithms cannot fully represent the actual status in the TME. Therefore, additional better-designed prospective studies should be performed to enhance the understanding in this field. Despite the promising results obtained in this study, experimental validation in both human tissue and cell lines should be performed in the future. However, the search for fully reliable biomarkers for HCC patient prognosis is still ongoing. We hope that our findings can provide specific insights into this field.

In conclusion, in this study, the immune microenvironment landscape of human hepatocellular carcinoma (HCC) and distinctive macrophage marker genes (MMGs) were comprehensively delineated with multiple single-cell and bulk RNA sequencing datasets. Also, a deeper analysis of macrophages based on machine learning defined three MMG-associated clusters and 12 hub genes. Furthermore, a novel scoring system named Macrosig score was constructed, characterized by different immune infiltration, mutation frequency and immune checkpoint expression, proved to have the potential to aid in efficacy prediction of immunotherapy and personalize of molecular therapeutic regimens.

HCC: Hepatocellular carcinoma;

LIHC: Liver hepatocellular carcinoma;

TME: Tumor microenvironment;

TAM: Tumor-associated macrophages;

CNV: Copy number variation;

TMB: Tumor mutation burden;

TCGA: The Cancer Genome Atlas;

GEO: Gene-Expression Omnibus;

DEGs: Differentially expressed genes;

UMAP:Uniform manifold approximation and projection;

PCA: Principal component analysis;

GSVA: Gene set variation analysis;

ssGSEA: Single-sample gene-set enrichment analysis;

ICIs: immune checkpoint inhibitors;

MMGs: Macrophage Marker Genes

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The scRNA-sequencing datasets generated in our study were acquired from TISCH (http://tisch.comp-genomics.org/). Publicly available gene expression data and complete clinical annotations from patients with LIHC were retrospectively obtained from the HCCDB (http://lifeome.net/database/hccdb/home.html), and the TCGA-LIHC data were obtained from the UCSC-XENA website (https://xenabrowser.net/datapages/). Further information is available from the corresponding author upon reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by Shanghai Rising-Star Program (Grant Nos. 23YF1441200), Fudan university Zhongshan Hospital Youth Fund (Grant Nos. 2021ZSQN49).

Author contributions

TL and XX equally contributed to this work. TL conceptualized and designed the study. TL, XX analyzed and interpreted the data, TL, XX and KN performed the bioinformatics analyses. TL, XX drafted the manuscript. ZR revised the manuscript for important intellectual content, LZ supervised the study and provided administrative and technical supports. All authors read and approved the final manuscript.

Acknowledgements

Not applicable.

Wen N, Cai Y, Li F, Ye H, Tang W, Song P, Cheng N: The clinical management of hepatocellular carcinoma worldwide: A concise review and comparison of current guidelines: 2022 update. Biosci Trends 2022, 16:20-30.
Parikh ND, Pillai A: Recent Advances in Hepatocellular Carcinoma Treatment. Clin Gastroenterol Hepatol 2021, 19:2020-2024.
Rizzo A, Ricci AD, Di Federico A, Frega G, Palloni A, Tavolari S, Brandi G: Predictive Biomarkers for Checkpoint Inhibitor-Based Immunotherapy in Hepatocellular Carcinoma: Where Do We Stand? Front Oncol 2021, 11:803133.
Villanueva A: Hepatocellular Carcinoma. N Engl J Med 2019, 380:1450-1462.
Gibney GT, Weiner LM, Atkins MB: Predictive biomarkers for checkpoint inhibitor-based immunotherapy. Lancet Oncol 2016, 17:e542-e551.
Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, Coussens LM, Gabrilovich DI, Ostrand-Rosenberg S, Hedrick CC, et al: Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med 2018, 24:541-550.
Ho DW, Tsui YM, Chan LK, Sze KM, Zhang X, Cheu JW, Chiu YT, Lee JM, Chan AC, Cheung ET, et al: Single-cell RNA sequencing shows the immunosuppressive landscape and tumor heterogeneity of HBV-associated hepatocellular carcinoma. Nat Commun 2021, 12:3684.
Oura K, Morishita A, Tani J, Masaki T: Tumor Immune Microenvironment and Immunosuppressive Therapy in Hepatocellular Carcinoma: A Review. Int J Mol Sci 2021, 22.
Noy R, Pollard JW: Tumor-associated macrophages: from mechanisms to therapy. Immunity 2014, 41:49-61.
Cassetta L, Pollard JW: Targeting macrophages: therapeutic approaches in cancer. Nat Rev Drug Discov 2018, 17:887-904.
Kuang DM, Zhao Q, Peng C, Xu J, Zhang JP, Wu C, Zheng L: Activated monocytes in peritumoral stroma of hepatocellular carcinoma foster immune privilege and disease progression through PD-L1. J Exp Med 2009, 206:1327-1337.
Qiu Y, Chen T, Hu R, Zhu R, Li C, Ruan Y, Xie X, Li Y: Next frontier in tumor immunotherapy: macrophage-mediated immune evasion. Biomark Res 2021, 9:72.
Nagaraju GP, Dariya B, Kasa P, Peela S, El-Rayes BF: Epigenetics in hepatocellular carcinoma. Semin Cancer Biol 2022, 86:622-632.
Chen H, Ye F, Guo G: Revolutionizing immunology with single-cell RNA sequencing. Cell Mol Immunol 2019, 16:242-249.
Rose TL, Weir WH, Mayhew GM, Shibata Y, Eulitt P, Uronis JM, Zhou M, Nielsen M, Smith AB, Woods M, et al: Fibroblast growth factor receptor 3 alterations and response to immune checkpoint inhibition in metastatic urothelial cancer: a real world experience. Br J Cancer 2021, 125:1251-1260.
Motzer RJ, Robbins PB, Powles T, Albiges L, Haanen JB, Larkin J, Mu XJ, Ching KA, Uemura M, Pal SK, et al: Avelumab plus axitinib versus sunitinib in advanced renal cell carcinoma: biomarker analysis of the phase 3 JAVELIN Renal 101 trial. Nat Med 2020, 26:1733-1741.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015, 43:e47.
Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD: The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012, 28:882-883.
Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, 3rd, Hao Y, Stoeckius M, Smibert P, Satija R: Comprehensive Integration of Single-Cell Data. Cell 2019, 177:1888-1902 e1821.
Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, Chak S, Naikawadi RP, Wolters PJ, Abate AR, et al: Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 2019, 20:163-172.
Wilkerson MD, Hayes DN: ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010, 26:1572-1573.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA, et al: Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013, 4:2612.
Hanahan D: Hallmarks of Cancer: New Dimensions. Cancer Discov 2022, 12:31-46.
Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012, 16:284-287.
Zhao P, Zhen H, Zhao H, Huang Y, Cao B: Identification of hub genes and potential molecular mechanisms related to radiotherapy sensitivity in rectal cancer based on multiple datasets. J Transl Med 2023, 21:176.
Zhang Q, He Y, Luo N, Patel SJ, Han Y, Gao R, Modak M, Carotta S, Haslinger C, Kind D, et al: Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma. Cell 2019, 179:829-845 e820.
Chen DS, Mellman I: Elements of cancer immunity and the cancer-immune set point. Nature 2017, 541:321-330.
Zhang X, Shi M, Chen T, Zhang B: Characterization of the Immune Cell Infiltration Landscape in Head and Neck Squamous Cell Carcinoma to Aid Immunotherapy. Mol Ther Nucleic Acids 2020, 22:298-309.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP: Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018, 28:1747-1756.
Markl F, Huynh D, Endres S, Kobold S: Utilizing chemokines in cancer immunotherapy. Trends Cancer 2022, 8:670-682.
Bule P, Aguiar SI, Aires-Da-Silva F, Dias JNR: Chemokine-Directed Tumor Microenvironment Modulation in Cancer Immunotherapy. Int J Mol Sci 2021, 22.
Propper DJ, Balkwill FR: Harnessing cytokines and chemokines for cancer therapy. Nat Rev Clin Oncol 2022, 19:237-253.
Curran MA, Montalvo W, Yagita H, Allison JP: PD-1 and CTLA-4 combination blockade expands infiltrating T cells and reduces regulatory T and myeloid cells within B16 melanoma tumors. Proc Natl Acad Sci U S A 2010, 107:4275-4280.
Llovet JM, Castet F, Heikenwalder M, Maini MK, Mazzaferro V, Pinato DJ, Pikarsky E, Zhu AX, Finn RS: Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol 2022, 19:151-172.
Xiong H, Mittman S, Rodriguez R, Moskalenko M, Pacheco-Sanchez P, Yang Y, Nickles D, Cubas R: Anti-PD-L1 Treatment Results in Functional Remodeling of the Macrophage Compartment. Cancer Res 2019, 79:1493-1506.
Deng L, He K, Pan Y, Wang H, Luo Y, Xia Q: The role of tumor-associated macrophages in primary hepatocellular carcinoma and its related targeting therapy. Int J Med Sci 2021, 18:2109-2116.
Chen J, Zheng DX, Yu XJ, Sun HW, Xu YT, Zhang YJ, Xu J: Macrophages induce CD47 upregulation via IL-6 and correlate with poor survival in hepatocellular carcinoma patients. Oncoimmunology 2019, 8:e1652540.
Xia Y, Rao L, Yao H, Wang Z, Ning P, Chen X: Engineering Macrophages for Cancer Immunotherapy and Drug Delivery. Adv Mater 2020, 32:e2002054.
Rizvi S, Wang J, El-Khoueiry AB: Liver Cancer Immunity. Hepatology 2021, 73 Suppl 1:86-103.
Yang C, Huang X, Liu Z, Qin W, Wang C: Metabolism-associated molecular classification of hepatocellular carcinoma. Mol Oncol 2020, 14:896-913.

No competing interests reported.

Additionalfile1.Fig.S1.tif
Additional file 1. Fig. S1: KM survival curve for each MMG.
Additionalfile2.Fig.S2.tif
Additional file 2. Fig. S2: Verification of the generally high expression of the 13 identified MMGs in macrophages in multiple liver cancer single-cell datasets from TISCH.
Additionalfile3.Fig.S3.tif
Additional file 3. Fig.S3: Comprehensive analysis of the 12 hub via the GSCA database. (A) The GSCA database was used to analyse the differential expression of the 12 hub genes in HCC via TCGA. The solid circles represent differential expression, and red represents high tumour expression. (B) Frequencies of deleterious mutations in HCC. (C) Oncoplot showing the mutation distribution of the top 10 mutated genes. (D) Pie chart summarizing the CNV profile of genes in the indicated tumour types.(E) Correlation between CNVs and gene expression. (F) Differences in methylation. (G) Correlation between methylation and mRNA expression for each gene.
Additionalfile4.Fig.S4.tif
Additional file 4. Fig. S4: Prediction of upstream miRNAs and transcription factors of the 12 hub genes via the RegNetwork database via Cytoscape.
Additionalfile5.Table.S1.xlsx
Additional file 5. Table. S1: List of all the differentially expressed marker genes.
Additionalfile6.Table.S2.xlsx
Additional file 6. Table. S2: MMG expression profiles in HCC in the univariate Cox regression analysis.
Additionalfile7.Table.S3.xlsx
Additional file 7. Table. S3: The Cox regression analysis with 761 DEGs.

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of a Macrophage Marker Gene Signature to Evaluate Immune Infiltration and Therapeutic Response in Hepatocellular Carcinoma Based on Single-cell and Bulk RNA Sequencing Data

Status:

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Materials and methods

Raw dataset acquisition

Bulk RNA-seq and scRNAseq data processing and analysis

Identification of MMG-associated clusters and estimation of immune cell infiltration

GSVA and functional enrichment analysis

Analyses of the hub genes and mutation landscapes of the MMG clusters

Construction and validation of a the prognosis-related Macrosig scoring system

Immunotherapy and targeted drug response prediction

Statistical analysis

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1