Comprehensive Analysis of Immune Prognostic‐Related Genes in the Tumor Microenvironment of Hepatocellular Carcinoma

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

Download PDF

Research article

Comprehensive Analysis of Immune Prognostic‐Related Genes in the Tumor Microenvironment of Hepatocellular Carcinoma

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

This work is licensed under a CC BY 4.0 License

Journal Publication

published 31 Mar, 2021

Read the published version in BMC Cancer →

You are reading this older preprint version

Read the latest preprint version →

Background: The percentage of death resulted from hepatocellular carcinoma (HCC) remains high worldwide, despite surgical and chemotherapy treatment. Immunotherapy offers great promise in the treatment of a rapidly expanding spectrum of HCC. Therefore, further exploration of the immune-related signatures in the tumor microenvironment, which plays a vital role in tumor initiation and progression for immunotherapy is currently needed.

Methods: In this present research, 866 immune-related difference expression genes (DEGs) were identified by integrating the DEGs between TCGA HCC and normal tissue and the immune genes from databases (Innate DB; Imm Port), and 144 candidate prognostic genes were defined through weighted gene co-expression network analysis (WGCNA).

Results: Seven prognostic immune-related DEGs were determined with LASSO Cox PH model, which was followed by constructing the ImmuneRiskScore model based on the prognostic immune-related DEGs. The prognostic index of the ImmuneRiskScore was validated then in the dependent dataset. Patients were divided into high- and low-risk groups according to ImmuneRiskScore. The difference in ImmuneRiskScore and infiltration of immune cells between groups was detected and the correlation analysis for immunotherapy biomarkers was further explored.

Conclusion:The ImmuneRiskScore of HCC could provide a prognostic signature and reflect immune characteristics within tumor microenvironment. Furthermore, it also may provide novel immunotherapy predictive biomarker for HCC patients in the near future.

Cancer Biology

Oncology

hepatocellular carcinoma

biomarker

prognostic

immune

Hepatocellular carcinoma (HCC) is one of the most common malignancies [1, 2]. With 5-year survival being 18%, HCC is the second most lethal tumor after pancreatic cancer[3] and the fourth leading cause of cancer-related mortality worldwide[4, 5]. Hepatocellular carcinoma is the major cancer type for primary liver cancers and its increase in deaths is a growing concern. However, general therapies such as radiotherapy and chemotherapy do not prolong overall survival (OS) significantly in HCC[6].

Thus, new strategies are required. Immunotherapy with checkpoint inhibitors is also emerging as an important treatment option. Immune checkpoint inhibitors (ICIs) is revolutionizing the clinical treatment landscape of multiple tumors, most notably advanced melanoma[7–13], non-small-cell lung cancer[14, 15] and renal cell carcinoma[16, 17]. Since HCC seems to be submissive to the programmed cell death protein 1 (PD1) pathway blockade[18], ICIs with approval likely for additional indications for HCC in the near future[6, 19, 20]. Though the remarkable progress that has been achieved, there is only a limited number of patients could benefit from ICIs[21]. Therefore, there is an urgent need for new, immune-based biomarkers distinguishing HCC patients more likely to have a better prognosis and further benefit from immunotherapy.

As an important element of the HCC tumor microenvironment, immune cells show clinicopathological significance in predicting prognosis and therapeutic efficacy[22–24]. It is an active field to investigate the characteristics of the tumor microenvironment functionally impact effect on immunotherapy.

In this study, we make full use of TCGA data and a priori immune-related genes to construct a prognostic immune risk score by WGCNA and Lasso cox model. We also analyzed the correlation for ImmuneRiskScore and different immune cells to elucidate mechanisms possible for the formation of the microenvironment. Last, we explored the correlation with other immune biomarkers and the potential to identify patients eligible for immunotherapy to improve the therapy effects. The flowchart as shown in Fig. 1.

Data download and procession

We downloaded the RNA sequencing expression profile(count and RPKM format) and clinical data of TCGA-HCC from the UCSC Xena data portal(https://xena.ucsc.edu/), which contains 50 normal samples and 374 tumor ones. The immune-related genes set contains 1052 immune genes downloaded from InnateDB (https://www.innatedb.ca/) (Tables S1) and 1811 immune-related genes downloaded from ImmPort[25](https://www.immport.org) (Tables S2). The expression profiling datasets GSE14520 as well as the clinical data were obtained by GEOquery[26] package of Bioconductor in R-3.5.2. GSE14520 contains 162 tumor samples respectively after removing the normal samples. Microarray probe ID mapped to gene symbols based on the GPL3921 platform (Affymetrix HT Human Genome U133A Array) and incorporated in the dataset matrix for each dataset. Eventually, the average of multiple probes computed that correspond to a single gene for each dataset individually employing in-house R scripts. The tumor mutation burden data was download from PanCancerAtla (https://www.cell.com/consortium/pancanceratlas).

Difference Expression Analysis

We identified differential expression genes (DEGs) based on limma[27] packages in R between Normal (n = 50) and cancer(n = 367) patients based on the raw counts of HCC gene expression data from TCGA. Empirical Bayes method was used in selecting significant DEGs based on "limma" package which was used in the standard comparison mode, and the threshold was set at p < 0.05 and |log2 fold change| >= 1.5.

Gene Ontology (GO) and pathway enrichment analysis of DEGs

As an ontology-based R package, clusterProfiler not only can automate the process of biological-term classification and the enrichment analysis of gene clusters, but also provides a visualization module for analysis results display[28]. In this present research, the clusterProfiler package was used in identifying and visualizing the GO terms (biological process, cellular component, and molecular function included) and KEGG pathways enriched by DEGs. We set P-value < 0.01 as the cut-off criterion, the significant adjustment method was BH and the cut-off criterion of q-value was also set 0.01.

WGCNA

We used fpkm format of TCGA datasets to construct WGCNA analysis[29]. First, Basic data preprocessing for handling missing data and removing outliers; Second, Choosing the soft-thresholding power by the pickSoftThreshold. It can calculate the scale-free topology fit index for several powers and provide the appropriate soft-thresholding power for network construction. Third, constructing a one-step automatic network and detecting modules. Turning adjacency into topological overlap to measure the network connectivity of a gene considered as the sum of its adjacency with all other genes for the generation of network. Hierarchical clustering function was used in classifying genes with a similar expression profile into modules[30]; Next, Selecting the key modules related to OS and OS time and visualized by Cytoscape[31]. In this present study, modules were chosen and were visualized. The correlation between MEs and clinical traits "survival" and "survival time" was calculated to identify the module related. And then, in the linear regression between gene expression and clinical information, we defined gene significance (GS) as the log10 transformation of the P-value (GS = lgP). In a module, module significance (MS) was defined as the average GS for all the genes. Hub genes were identified with high clinical traits significance (> 0.1) and high intramodular connectivity (> 0.5) in interesting modules and were selected as the candidates to be further analysed and validated.

Gene set enrichment analysis of hub genes

We used gprofiler2 (https://CRAN.R-project.org/package=gprofiler2) to perform over-representation analysis on input HCC hub gene list. Based on high-quality up-to-date data across different evidence types, g:Profiler can provide a reliable service[32]. It maps these immune genes to known functional information sources and detects statistically significantly enriched terms. We included pathways from KEGG(https://www.genome.jp/kegg/), Reactome (https://reactome.org/), and CORUM (http://mips.helmholtz-muenchen.de/corum/). This method was done with the hypergeometric test followed by FDR correction for multiple testing.

Construction and Validation of an Immunoscore Prognostic Model

With the univariable Cox proportional hazards regression model in 'survival' package, we made a calculation of the hazard proportions for DEGs of the HCC cohort. DEGs with significance p value < 0.05 were analyzed, and we re-computed these DEGs with survival risk for Cox regression models by glmnet package[32]to select several more important prognostic genes among DEGs. As an R package, glmnet fits a generalized linear model via penalized maximum likelihood. The regularization path is computed for the lasso by setting the regularization parameter lambda as 1. To predict patient survival, a formula for the ImmuneRiskScore model was established as follows:

Immune Risk Score =\(\sum \left(Normalized expression value of gene Gi\text{ × Lasso cox coefficient of gene Gi}\right)\)

Estimation of Immune Cell Type Fractions and estimate score

CIBERSORT[33] can characterize cell composition of complex tissues from their gene expression profiles. In this study, we used CIBERSORT and LM22 reference gene expression matrix to quantify cell composition in diverse HCC samples. Analysis of normalized gene expression data were carried out by using the CIBERSORT algorithm, running with 1,000 permutations. Calculation of immune and stromal scores were performed with ESTIMATE, which was an algorithm could provid scores for the level of stromal cells present, and the infiltration level of immune cells in tumor tissues[34].

Statistical Analysis

The survival curve for the hub immune genes was created by the Kaplan-Meier method and the statistical significance of difference was judged by the Log-rank test. The receiver operating characteristic (ROC) curve was applied to describe the sensitivity and specificity of survival prediction based on the Immune Risk Score, and the pROC package was utilized to quantify the area under the curve (AUC). Nonparametric Mann-Whitney-Wilcoxon Test was used to compare the data from different groups and Pearson's chi-square test was performed in measuring the level of significance for association amongst variables. All statistical analyses were conducted using R software. All of the reported P values were two-tailed, and p < 0.05 was considered to be statistically significant.

Identification of DEGs related to immune

From the TCGA database, we obtained the expression profiles of 417 HCC samples, which contained 367 tumor samples and 50 nontumor ones after data preprocessing. Totally 7194 genes were identified as DEGs with the threshold of P < 0.05 and |log2FC| > 1.5, containing 3657 genes upregulated and 3537 genes downregulated (Fig. 2A, Tables S3). The samples could be well clustered for normal and tumor groups when the top 200 DEGs were selected for unsupervised hierarchical clustering (Fig. 2B). To obtained the immune-related DEGs for HCC samples, we selected the genes both be immune-related genes and difference expressed for HCC. We total collected 2542 human immune-related genes composed by the 1811 immune-related genes from Immport and 1052 innate immune-related genes from InnatedDB. Finally, there were 866 genes within all the DEGs when intersecting with immune genes (Fig. 1C, Tables S4), and These immune-related DEGs were then chosen for further analysis.

The further GO enrichment analysis of these 866 immune-related DEGs showed that 1008 go terms including 891 biological processes (BPs), 39 molecule functions (MFs), and 78 cellular components (CCs) were significantly enriched in (Tables S5). As shown in Fig. 3A, the top 10 go terms were showed by classification. Leukocyte migration, positive regulation of cytokine production, cell chemotaxis, leukocyte cell to cell adhesion, and T lymphocyte activation are the top enriched biological process; the enriched cellular components including collagen − containing extracellular matrix, external side of plasma membrane, collagen trimer, secretory granule membrane, extracellular matrix component; the enriched molecule functions contains structural constituent of extracellular matrix, carbohydrate binding, receptor ligand activity, cytokine activity, glycosaminoglycan binding. More details about the top GO enrichment analysis results were showed by overview their genes in Fig. 3B, C, D. These circle plots significantly could help us to understand the function of DEGs in these enriched terms. And most of the time, it is more meaningful to further represent genes with various functions. Furthermore, analysis of KEGG pathway enrichment (Tables S6) indicated that significantly enriched pathways include interaction between cytokine and cytokine receptor, Viral protein interaction with cytokine and cytokine receptor, Cell adhesion molecules (CAMs), Chemokine signaling pathway, Malaria (Fig. 3E). The top 10 pathways were showed and interact with their assigned genes in Fig. 3F. The KEGG pathways analysis results were different from GO enrichment results, indicating that the immune microenvironment of HCC was sophisticated.

Weighted Co-expression Network Construction And Key Modules Identification

The 866 immune-related DEGs and 367 HCC tumor samples were used in constructing gene co-expression network. After checked the missing values, we detected the outlier samples by hierarchical clustering (Figure S1), and the dendrogram shows 5 outlier samples that will be removed from the analysis. In this study, β = 4 (scale-free R2 = 0.85) was selected to be the soft-thresholding parameter to ensure a scale-free network. As shown in Fig. 4A, we performed network topology analysis for thresholding powers from 1 to 20, 4 was the lowest power for the scale-free topology fit index on 0.85. We obtained gene clustering tree by using hierarchical clustering of TOM-based dissimilarity and identified 6 modules (Fig. 4B, Table 1). To select the clinically significant modules, we used WGCNA to calculate the correlations between the external clinical information and gene modules. As shown in Fig. 4C, the Green module was found to be the highest one associated with overall survival, and the green and brown eigengenes are highly related (Figure S2).

Table 1

The number of genes contained by modules
Module	The number of genes
Blue	302
Brown	147
Green	139
Red	86
Turquoise	389
Yellow	139
Grey	254

We visualized the green module as networks by Cytoscape and selected the top 100 gene pairs by sorting the weight of gene pairs. As shown in Fig. 4D, it indicates that some genes, such as PDLIM7, EHHADH, DMGDH, and CYP8B1 with larger size have higher node degree. Finally, we screened 144 immune-hub genes with high gene significance with OS and OS time (> 0.1), and high relationship with interesting green module (> 0.5) (Table S7).

Statistical gene set enrichment analysis about the 144 immune-hub genes was performed to find over-representation of functions from biological pathways like KEGG, Reactome, and complexes in CORUM, etc. This is done with the hypergeometric test followed by correction for multiple testing. The results showed that immune-hub genes were enriched significantly in 52 pathways and complexes (p < 0.05) (Fig. 4E; Tables S8), including KEGG pathways such as TNF signaling pathway, Metabolic pathways, Reactome pathways such as Phenylalanine and tyrosine metabolism, and CORUM complexes such as PLAUR − PLAU complex, IGF2R − PLAUR − PLAU complex, MAK − ACTR − AR complex, IGF2R − PLG − PLAU − PLAUR − LTGFbeta1 complex. These findings show that these hub genes not only affect the metabolic, apoptosis, and cell survival as well as inflammation and immunity of HCC, but also plays a pivotal role in the protein complexes of immune cells.

Establishment of the lasso cox-based prognostic gene signature

We subsequently revealed that 108 of the 144 immune-hub genes were significantly related to OS through univariate Cox regression analysis (results in Tables S9). Then, we performed lasso-penalized Cox analysis to further narrow the hub genes (Fig. 5AB). Seven genes were identified and utilized thereafter in constructing an ImmuneRiskScore model to evaluate the prognostic value of CRC patients. The seven genes identified and their cox coefficient were showed in Table 2 including Secreted Phosphoprotein 2(SPP2), Glucose-6-Phosphatase Catalytic Subunit(G6PC), Phospholipase B Domain Containing 1(PLBD1), ETS Variant Transcription Factor 4(ETV4), Phosphofructokinase Platelet(PFKP), G Protein Subunit Alpha Z(GNAZ), Asparaginase And Isoaspartyl Peptidase 1(ASRGL1). The GO enrichment analysis shows that these seven genes were enrichment in several molecular functions such as asparaginase activity, beta-aspartyl-peptidase activity, 6-phosphofructokinase activity, glucose-6-phosphatase activity, and sugar-terminal-phosphatase activity. The formula for the ImmuneRiskScore model was described in the part of Materials and Methods.

Next, we divided HCC patients into high- and low-score groups, based on the lasso cox hub genes and ImmuneRiskScore, according to the optimal cut-off got from survminer package. The results indicated that five genes (PLBD1, ETV4, PFKP, GNAZ, ASRGL1) are risk factors, and that high-score samples had a worse OS than those with low score and the SPP2 and G6PC are protective factors (Fig. 5C). The last figure in Fig. 5C showed the prognostic accuracy of ImmuneRiskScore (95% CI for HR: 0.48 (0.33–0.68), log-rank test p < 0.0001). Additionally, the result of multivariate Cox regression analyses showed that the predictive value of the ImmuneRiskScore for HCC patients is not associated with common clinical variables (Tables S10).

Table 2

The result of Lasso regression
Genes	Coef.
SPP2	-0.27112
G6PC	-0.12328
PLBD1	0.66909
ETV4	0.359013
PFKP	0.402211
GNAZ	0.126409
ASRGL1	0.795796

Validation Of The Immuneriskscore In Tcga Crc Cohort

To further investigate the prognostic value of ImmuneRiskScore, we conducted a validation analysis in another GEO cohort (GSE14520). The dataset was categorized into two groups based on ImmuneRiskScore; the results of analysis indicate that there is a significant prognostic value between ImmuneRiskScore and OS as well as recurrence, the high-score group had a worse OS rate than the low-score one, not merely in the TCGA cohort (Fig. 6AB). Figure 6C showed the prognostic accuracy of ImmuneRiskScore in the GEO dataset, which was considered as a continuous variable during investigation. The area under the ROC curves (AUC) of the prognostic model for OS was 0.608 at 1 year, 0.614 at 3 years, 0.620 at 5 years. These results indicated that ImmuneRiskScore is a good model to predict survival.

Stromal and Immune cell infiltration landscapes between high and low ImmuneRiskScore

By using ESTIMATE algorithm, we estimated the infiltrating cells and tumor purity of tumor tissue. Stromal score represented the presence of stromal cells in tumor tissue, and immune score indicated the infiltration of immune cells in tumor tissue, and combination of them stood for a measurement of tumor purity (Tables S11). We then distinguished differences in stromal and immune scores between high- and low-risk patients with LUSC. As shown in Fig. 7(A, B), the immune and stromal score in high-level patients were both significantly (Wilcox p < 0.05) higher than low-level HCC patients indicate the overall presence of immune and stromal levels are associated with Immune Risk Score.

We further estimated the ratio of 22 immune cell categories in HCC patients using the CIBERSORT method, the results are shown in Tables S12. The proportion distribution of immune cells is heterogeneity not only between the high and low-levels samples of ImmuneRiskScore but among the HCC samples (Figure S3). We made a comparison of the relative proportions of 22 immune cells between low and high ImmuneRiskScore in HCC patients, and found significant differences in B cells memory, Plasma cells, T cells CD4 memory activated, T cells regulatory Tregs, NK cells resting, Monocytes, Macrophages M0, Dendritic cells resting, Mast cells resting, Neutrophils(Fig. 7C). Furthermore, we calculated the absolute immune infiltrated score of these 22 immune cells by combined with tumor purity. A similar comparison of the relative proportions of 22 types of immune cells between low and high ImmuneRiskScore was conducted and significant differences were found in B cells memory, Plasma cells, T cells CD4 memory activated, T cells follicular helper, T cells regulatory Tregs, NK cells resting, NK cells activated, Macrophages M0, Macrophages M2, Dendritic cells resting, Mast cells resting, Neutrophils. It can be easily obtained that several immune cells were associated with ImmuneRiskScore both in relative and absolute types including B cells memory, Plasma cells, T cells CD4 memory activated, T cells regulatory Tregs, NK cells resting, Macrophages M0, Dendritic cells resting, Mast cells resting, Neutrophils. These results also demonstrated that the variations in levels of immune cells in tumor environment might have a correlationship with the overall survival in HCC samples.

Correlations between the ImmuneRiskScore and immune biomarkers

The discovery of broad immune biomarkers of the TME could effectively predict clinical benefit to ICIs. We next introduced a few important immune biomarkers including PDL1, PD1, PDL2, CTLA4, CYT, and IFN gamma. Among these biomarkers, the immune checkpoint genes including PDL1, PD1, PDL2, CTLA4 are worth attentional for Immune checkpoint inhibitors which are revolutionizing the clinical treatment landscape. Cytolytic activity (CYT) focuses on cytotoxic T cells (CTL) and natural killer cells (NK) for their powerful ability to lyse tumor cells[35]. CYT is measured based on the geometric mean of expression of granzyme A (GZMA) and perforin (PRF1). Interferon-γ(IFN-gamma) signature is a key cytokine to activate the PD-1 signaling axis by directly upregulating the ligands PD-L1 and PD-L2 mainly in tumor cell, which was produced by T cells activated, NK and NKT cells[36, 37]. Tumor mutational burden (TMB) is defined as the number of nonsynonymous mutations detected in each megabase sequenced. TMB is proved to be associated with improved responses to checkpoint blockade, in some tumors such as melanoma[38] and non-small-cell lung cancer[39, 40]. An experimentally determined pan-fibroblast TGF-β response signature (Pan-F-TBRS) showed elevated mean expression in non-responders and decreased overall survival particularly in patients with mUC[41]. transforming growth factor β (TGF-β) signalling in fibroblasts is documented as a pleiotropic cytokine having relationship with poor prognosis in multiple tumor categories[42, 43], and it is thought to be vital in advanced cancers in promotion of immunosuppression, angiogenesis, metastasis, tumor cell epithelial to mesenchymal transition (EMT), fibroblast activation and desmoplasia[44–46].

We further explored the relationship between ImmuneRiskScore and these immune biomarkers (Tables S13, Tables S14) by person which satisfying a bivariate normal distribution (Pearson result in Tables S15). As shown in Fig. 8, ImmuneRiskScore were significantly positive correlated with both several immune inflammed biomarkers(PDL1: r = 0.31;p = 1.63e-09,95% CI: 0.21–0.40; PD-1: r = 0.35;p = 8.45e-12,95% CI: 0.25–0.43; PD-L2: r = 0.27;p = 1.02e-07,95% CI: 0.18–0.37; CTLA-4: r = 0.42;p = 6.78e-05,95% CI: 0.33–0.50; CYT: r = 0.21;p = 5.41e-17,95% CI: 0.11–0.30, IFN gamma: r = 0.42;p = 0.00037,95% CI: 0.17–0.36) and Pan-F-TBRS (r = 0.27;p = 2.51e-07,95% CI: 0.33–0.50). All of the significant p value for correlation was smaller than 0.001, suggesting that ImmuneRiskScore might be a potential biomarker for immunotherapy especially for Immune checkpoint inhibitors. In addition, there is no correlation between ImmuneRiskScore and TMB, indicating TMB is an independent factor mediating the TME in these two groups.

Increasing evidence indicates that immune-related biomarker is associated with prognostic for various cancer types[47–49]. Moreover, prognostic biomarkers effectively guiding cancer therapy especially immunotherapy are still in need. Therefore, we constructed ImmuneRiskScore that contributes to overall survival by investigating the tumor microenvironment of HCC. Based on the gene expression data of HCC from TCGA, we identified a prognostic signature of seven immune hub genes (SPP2, G6PC, PLBD1, ETV4, PFKP, GNAZ, ASRGL1) through combining WGCNA, univariate Cox regression analysis, and LASSO PH model. Then, we constructed a seven-gene risk scoring system and classified HCC patients into two risk groups with significantly different survival rates. Successfully validation of the prognostic performance of the risk scoring model was conducted in an independent set from GEO. It shows that the seven-genes ImmuneRiskScore are promising prognostic biomarkers for HCC and may function importantly in the TME of HCC.

Generally, the ImmuneRiskScore is mainly the result of the presence of TME which including the complex interaction between cancer cells, stromal cells, and immune cells. Under this condition, we aimed to further explore the relationship between them. We found that the relative or absolute infiltration levels of B cells memory, Plasma cells, T cells CD4 memory activated, T cells regulatory Tregs, NK cells resting, Macrophages M0, Dendritic cells resting, Mast cells resting, Neutrophils are significantly correlated with ImmuneRiskScore, indicated that these cells are likely to have prognostic values. Such as infiltration of macrophages in solid tumors is associated with poor prognosis and correlates with chemotherapy resistance in most cancers[50]. It is noteworthy that these cell types didn’t conclude most of T cell compartment, which the clinical response strategies have largely focused on. In fact, other immune cells may also contribute to anti-tumour immunity[51–53]. Especially memory B cells also has potential role in the response to ICB treatment[54].

With an aim to effectively predict clinical benefit to checkpoint inhibitor strategies, we performed a wider exploration of active innate and adaptive immune responses within the tumor microenvironment by gene expression profiling. The predictive value of our immune relation score is envisaged for the positive correlation with PDL1, PD-1, PD-L2, CTLA-4, CYT, IFN gamma and Pan-F-TBRS. These positive related biomarkers are involved in the proinflammatory cytokines related inflammation microenvironment of tumor [55, 56] and the TGF-β signal pathway related excluded microenvironment of tumor[41]. Corresponds to the previous infiltration results, the inflammation of TME can be measured by the cellular content of the tumor, for example, infiltrating T and B cells. Inflamed tumors also contain a broad proinflammatory cytokines profile and a type I interferon signature indicating activation of innate immune response. While TGF-β can drive an excluded phenotype in TME as its impact on stromal cells, and prevent T cell penetration into the centre of the tumor, [41]. These results indicate that anti-tumor immunological effect is a bidirectional and dynamic system in the tumor microenvironment (TME). This biomarker analysis should help us to unravel the complexities of the interaction and molecular mechanisms between cancer and the host immune system.

In general, we gained a comprehensive insight into the TME of HCC and created a prognostic immune-related score that might become potential prognostic and predictive biomarkers.

AbbreviationsGEO: Gene Expression Omnibus; TCGA: The Cancer Genome Atlas; HCC: Hepatocellular carcinoma; TME: Tumor microenvironment; TMB: Tumor mutation burden; DEGs: difference expression genes; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: gene ontology; PDL1: Programmed cell death receptor ligand 1; PD-1: Programmed cell death receptor 1; PD-L2: Programmed cell death receptor ligand 1; CTLA-4: Cytotoxic lymphocyte antigen 4; CYT: Cytolytic Activity; IFN: Interferon; TGF-β: Transforming growth factor β; ICIs: Immune checkpoint inhibitors; ICB: immune checkpoint blockade; WGCNA: weighted gene co-expression network analysis; LASSO Cox PH model: L1-penalized least absolute shrinkage and selection operator (LASSO) Cox proportional hazards (PH) model; OS: Overall survival

Figure titles and legends

Figure 1 Flowchart

Figure 2 Identification of DEGs in patients with HCC for normal and tumor samples. (A) volcano plot of 7194 DEGs, upper-left and upper-right blue dots stand for genes down- and up-regulated in HCC, respectively. (B) Unsupervised hierarchical clustering heatmap for top 200 DEGs ranked by foldchange. Firebrick: up-regulated DEGs; navy blue: down-regulated DEGs. (C) Venn plot for Immune gene set and DEGs.

Figure 3 GO enrichment analysis result of immune-related DEGs. The bubble plot for GO enrichment analysis (A) and KEGG (B). The top 10 terms by adjusted p-values were selected to show. The x-axis of the plot represents the z-score. The -log (adjusted p-value) is displayed on the y-axis which corresponding to the significance of the term. The area of the plotted circles is proportional to the number of genes assigned to the term. Each circle is colored according to its category and labeled with the ID. The yellow lines represent a threshold for the displayed labels (adjusted p-value < 0.05).

Figure 4 WGCNA and enrichment analysis. A: The scale-free fit index analysis for various soft-thresholding powers, cut-off > 0.85. The mean connectivity analysis for various soft-thresholding powers. In all, 4 was the fittest power value (β). B: The cluster dendrogram of DEGs. In the figure, each branch represents one gene, and every color below indicates one co-expression module. C: Heatmap of the correlation between module eigengenes and the clinical traits of HCC. The green module was correlated mostly with OS status and time. D: The top 50 hub genes in the green module. Genes are represented by nodes, and node size is correlated with connectivity of the gene by degree, and edge size is correlated with weight. E: Visualization and table of the enrichment results. left: The x-axis stands for the functional terms that are grouped and color-coded according to data sources. The y-axis indicates the adjusted p-values in negative log10 scale. Every circle is one term and is sized according to the term size. right: the table of interesting enrichment results.

Figure 5 LASSO Cox regression analysis process for HCC immune-related genes. A: The solution path plot of each independent variable. Lateral axis and longitudinal axis stand for the lambda value and independent variable coefficient, respectively. Each curve corresponds to a variable. B: The confidence interval under each lambda. It contains the cross-validation curve (red dotted line), and upper and lower standard deviation curves along the λ sequence (error bars). Two selected λ’s are shown by the vertical dotted lines. LASSO, least absolute shrinkage and selection operator. C: Prognosis analysis of the selected Lasso cox hub genes and ImmuneRiskScore.

Figure 6 Validation of ImmuneRiskScore model; Patients were grouped by low or high ImmuneRiskScore (low- and high-score). Kaplan-Meier curves for OS (A) and recurrence (B), time-dependent ROC curve analysis of the ImmuneRiskScore (C) in GSE14520 cohort.

Figure 7 Landscape of immune microenvironment between low and high ImmuneRiskScore groups. The landscape of ImmuneScore (A) and StromalScore (B) of HCC patients in low and high ImmuneRiskScore groups. (C) The relative ratio of immune cell categories in low and high ImmuneRiskScore groups. (D) The absolute proportion of immune cell categories in low and high ImmuneRiskScore groups. Comparison of the statistical differences between the two groups were performed through Wilcox rank-sum test. Each boxplot is labeled with asterisks for p values. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001)

Figure 8 Correlation scatterplots between ImmuneRiskScore and immune check block (ICB) associated biomarkers, combined with density plot of expression distribution. The ICB biomarkers including PDL1, PD1, PDL2, CTLA4, CYT, IFN gamma, Pan_F_TBRS and TMB.

GEO: Gene Expression Omnibus; TCGA: The Cancer Genome Atlas; HCC: Hepatocellular carcinoma; TME: Tumor microenvironment; TMB: Tumor mutation burden; DEGs: difference expression genes; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: gene ontology; PDL1: Programmed cell death receptor ligand 1; PD-1: Programmed cell death receptor 1; PD-L2: Programmed cell death receptor ligand 1; CTLA-4: Cytotoxic lymphocyte antigen 4; CYT: Cytolytic Activity; IFN: Interferon; TGF-β: Transforming growth factor β; ICIs: Immune checkpoint inhibitors; ICB: immune checkpoint blockade; WGCNA: weighted gene co-expression network analysis; LASSO Cox PH model: L1-penalized least absolute shrinkage and selection operator (LASSO) Cox proportional hazards (PH) model; OS: Overall survival

Funding

Not applicable.

Availability of data and materials

The datasets generated and analyzed during the current study are available in the TCGA(https://xena.ucsc.edu/) and GEO repository with the dataset number GSE14520 , (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14520).

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

Peiming Guo and Guangjun Shi participated in the conception and design of the study. Weike Gao and Luan Li partially participated in the design of the study, final preparation and revision of the paper. Shuai Hong and Dongliang Wang performed bioinformatics analyses. Chengzhen Li, Guanying Yu, Lei Zhang, Dongsheng Zhang and Caiyun Liu drafted the paper. All authors read and approved this final manuscript.

Acknowledgments

Not applicable.

Siegel RL, Miller KD, Jemal A. Cancer statistics, 2016. CA Cancer J Clin. 2016;66(1):7–30.
El-Serag HB. Hepatocellular carcinoma. N Engl J Med. 2011;365(12):1118–27.
Annual Report to. the Nation on the Status of Cancer. Jnci Journal of the National Cancer Institute 2008.
Higginson J. International Agency for Research on Cancer. Encyclopedia of Toxicology. 1969;22(12):517–22.
WHO Press
McGuire S. World Cancer Report 2014. Geneva, Switzerland: World Health Organization, International Agency for Research on Cancer, WHO Press, 2015. Advances in Nutrition An International Review Journal 2016, 7(2):418–419.
Greten TF, Lai CW, Li G, Staveley-O'Carroll KF. Targeted and Immune-Based Therapies for Hepatocellular Carcinoma. Gastroenterology. 2019;156(2):510–24.
Topalian SL, Sznol M, McDermott DF, Kluger HM, Carvajal RD, Sharfman WH, Brahmer JR, Lawrence DP, Atkins MB, Powderly JD, et al. Survival, durable tumor remission, and long-term safety in patients with advanced melanoma receiving nivolumab. J Clin Oncol. 2014;32(10):1020–30.
Hamid O, Robert C, Daud A, Hodi FS, Hwu WJ, Kefford R, Wolchok JD, Hersey P, Joseph RW, Weber JS, et al. Safety and tumor responses with lambrolizumab (anti-PD-1) in melanoma. N Engl J Med. 2013;369(2):134–44.
Wolchok JD, Kluger H, Callahan MK, Postow MA, Rizvi NA, Lesokhin AM, Segal NH, Ariyan CE, Gordon RA, Reed K, et al. Nivolumab plus ipilimumab in advanced melanoma. N Engl J Med. 2013;369(2):122–33.
Robert C, Thomas L, Bondarenko I, O'Day S, Weber J, Garbe C, Lebbe C, Baurain JF, Testori A, Grob JJ, et al. Ipilimumab plus dacarbazine for previously untreated metastatic melanoma. N Engl J Med. 2011;364(26):2517–26.
Robert C, Schachter J, Long GV, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil C, Lotem M, et al. Pembrolizumab versus Ipilimumab in Advanced Melanoma. N Engl J Med. 2015;372(26):2521–32.
Robert C, Long GV, Brady B, Dutriaux C, Maio M, Mortier L, Hassel JC, Rutkowski P, McNeil C, Kalinka-Warzocha E, et al. Nivolumab in previously untreated melanoma without BRAF mutation. N Engl J Med. 2015;372(4):320–30.
Larkin J, Chiarion-Sileni V, Gonzalez R, Grob JJ, Cowey CL, Lao CD, Schadendorf D, Dummer R, Smylie M, Rutkowski P, et al. Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. N Engl J Med. 2015;373(1):23–34.
Muller M, Schouten RD, De Gooijer CJ, Baas P. Pembrolizumab for the treatment of non-small cell lung cancer. Expert Rev Anticancer Ther. 2017;17(5):399–409.
Brahmer J, Reckamp KL, Baas P, Crino L, Eberhardt WE, Poddubskaya E, Antonia S, Pluzanski A, Vokes EE, Holgado E, et al. Nivolumab versus Docetaxel in Advanced Squamous-Cell Non-Small-Cell Lung Cancer. N Engl J Med. 2015;373(2):123–35.
Motzer RJ, Escudier B, McDermott DF, George S, Hammers HJ, Srinivas S, Tykodi SS, Sosman JA, Procopio G, Plimack ER, et al. Nivolumab versus Everolimus in Advanced Renal-Cell Carcinoma. N Engl J Med. 2015;373(19):1803–13.
Motzer RJ, Rini BI, McDermott DF, Redman BG, Kuzel TM, Harrison MR, Vaishampayan UN, Drabkin HA, George S, Logan TF, et al. Nivolumab for Metastatic Renal Cell Carcinoma: Results of a Randomized Phase II Trial. J Clin Oncol. 2015;33(13):1430–7.
Sia D, Jiao Y, Martinez-Quetglas I, Kuchuk O, Villacorta-Martin C, Castro de Moura M, Putra J, Camprecios G, Bassaganyas L, Akers N, et al. Identification of an Immune-specific Class of Hepatocellular Carcinoma, Based on Molecular Features. Gastroenterology. 2017;153(3):812–26.
Makarova-Rusher OV, Medina-Echeverz J, Duffy AG, Greten TF. The yin and yang of evasion and immune activation in HCC. J Hepatol. 2015;62(6):1420–9.
Breous E, Thimme R. Potential of immunotherapy for hepatocellular carcinoma. J Hepatol. 2011;54(4):830–4.
El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, Kim TY, Choo SP, Trojan J, Welling THR, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet. 2017;389(10088):2492–502.
Godfrey DI, Le Nours J, Andrews DM, Uldrich AP, Rossjohn J. Unconventional T Cell Targets for Cancer Immunotherapy. Immunity. 2018;48(3):453–73.
Saito T, Nishikawa H, Wada H, Nagano Y, Sugiyama D, Atarashi K, Maeda Y, Hamaguchi M, Ohkura N, Sato E, et al. Two FOXP3(+)CD4(+) T cell subpopulations distinctly control the prognosis of colorectal cancers. Nat Med. 2016;22(6):679–84.
Khemlina G, Ikeda S, Kurzrock R. The biology of Hepatocellular carcinoma: implications for genomic and immune therapies. Mol Cancer. 2017;16(1):149.
Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, Hu Z, Zalocusky KA, Shankar RD, Shen-Orr SS, et al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data. 2018;5:180015.
Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007;23(14):1846–7.
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 research. 2015;43(7):e47.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: a journal of integrative biology. 2012;16(5):284–7.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.
A general framework for weighted gene
Zhang B, Horvath S. A general framework for weighted gene.
co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:Article17.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, Vilo J. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic acids research. 2019;47(W1):W191–8.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
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.
Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1–2):48–61.
Ikeda H, Old LJ, Schreiber RD. The roles of IFN gamma in protection against tumor development and cancer immunoediting. Cytokine Growth Factor Rev. 2002;13(2):95–109.
Fuchs CS, Doi T, Jang RW, Muro K, Satoh T, Machado M, Sun W, Jalal SI, Shah MA, Metges JP, et al. Safety and Efficacy of Pembrolizumab Monotherapy in Patients With Previously Treated Advanced Gastric and Gastroesophageal Junction Cancer: Phase 2 Clinical KEYNOTE-059 Trial. JAMA Oncol. 2018;4(5):e180013.
Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, Sucker A, Hillen U, Foppen MHG, Goldinger SM, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015;350(6257):207–11.
Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, Lee W, Yuan J, Wong P, Ho TS, et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015;348(6230):124–8.
Hellmann MD, Callahan MK, Awad MM, Calvo E, Ascierto PA, Atmaca A, Rizvi NA, Hirsch FR, Selvaggi G, Szustakowski JD, et al. Tumor Mutational Burden and Efficacy of Nivolumab Monotherapy and in Combination with Ipilimumab in Small-Cell Lung Cancer. Cancer Cell. 2018;33(5):853–61 e854.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8.
Lin RL, Zhao LJ. Mechanistic basis and clinical relevance of the role of transforming growth factor-beta in cancer. Cancer Biol Med. 2015;12(4):385–93.
Calon A, Lonardo E, Berenguer-Llergo A, Espinet E, Hernando-Momblona X, Iglesias M, Sevillano M, Palomo-Ponce S, Tauriello DV, Byrom D, et al. Stromal gene expression defines poor-prognosis subtypes in colorectal cancer. Nat Genet. 2015;47(4):320–9.
Massague J. TGFbeta in Cancer. Cell. 2008;134(2):215–30.
Derynck R, Zhang YE. Smad-dependent and Smad-independent pathways in TGF-beta family signalling. Nature. 2003;425(6958):577–84.
Flavell RA, Sanjabi S, Wrzesinski SH, Licona-Limon P. The polarization of immune cells in the tumour environment by TGFbeta. Nat Rev Immunol. 2010;10(8):554–67.
Shen S, Wang G, Zhang R, Zhao Y, Yu H, Wei Y, Chen F. Development and validation of an immune gene-set based Prognostic signature in ovarian cancer. EBioMedicine. 2019;40:318–26.
Choi H, Na KJ. Integrative analysis of imaging and transcriptomic data of the immune landscape associated with tumor metabolism in lung adenocarcinoma: Clinical and prognostic implications. Theranostics. 2018;8(7):1956–65.
Teschendorff AE, Miremadi A, Pinder SE, Ellis IO, Caldas C. An immune response gene expression module identifies a good prognosis subtype in estrogen receptor negative breast cancer. Genome Biol. 2007;8(8):R157.
Cassetta L, Pollard JW. Targeting macrophages: therapeutic approaches in cancer. Nat Rev Drug Discov. 2018;17(12):887–904.
Sarvaria A, Madrigal JA, Saudemont A. B cell regulation in cancer and anti-tumor immunity. Cell Mol Immunol. 2017;14(8):662–74.
Fridman WH, Zitvogel L, Sautes-Fridman C, Kroemer G. The immune contexture in cancer prognosis and treatment. Nat Rev Clin Oncol. 2017;14(12):717–34.
Tsou P, Katayama H, Ostrin EJ, Hanash SM. The Emerging Role of B Cells in Tumor Immunity. Cancer Res. 2016;76(19):5597–601.
Helmink BA, Reddy SM, Gao J, Zhang S, Basar R, Thakur R, Yizhak K, Sade-Feldman M, Blando J, Han G, et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature. 2020;577(7791):549–55.
Taube JM, Anders RA, Young GD, Xu H, Sharma R, McMiller TL, Chen S, Klein AP, Pardoll DM, Topalian SL, et al. Colocalization of inflammatory response with B7-h1 expression in human melanocytic lesions supports an adaptive resistance mechanism of immune escape. Sci Transl Med. 2012;4(127):127ra137.
Gajewski TF, Schreiber H, Fu YX. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol. 2013;14(10):1014–22.

Download PDF

Journal Publication

published 31 Mar, 2021

Read the published version in BMC Cancer →

Editorial decision: Major revision
17 Aug, 2020
Review #4 received at journal
15 Aug, 2020
Review #3 received at journal
15 Aug, 2020
Review #5 received at journal
12 Aug, 2020
Review #2 received at journal
08 Aug, 2020
Reviewer #5 agreed at journal
22 Jul, 2020
Review #1 received at journal
22 Jul, 2020
Reviewer #4 agreed at journal
21 Jul, 2020
Reviewer #3 agreed at journal
19 Jul, 2020
Reviewers invited by journal
17 Jul, 2020
Reviewer #1 agreed at journal
17 Jul, 2020
Reviewer #2 agreed at journal
17 Jul, 2020
Editor assigned by journal
14 Jul, 2020
First submitted to journal
13 Jul, 2020
Submission checks completed at journal
13 Jul, 2020
Editor invited by journal
13 Jul, 2020

You are reading this older preprint version

Read the latest preprint version →

Comprehensive Analysis of Immune Prognostic‐Related Genes in the Tumor Microenvironment of Hepatocellular Carcinoma

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Materials And Methods

Data download and procession

Difference Expression Analysis

Gene Ontology (GO) and pathway enrichment analysis of DEGs

WGCNA

Gene set enrichment analysis of hub genes

Construction and Validation of an Immunoscore Prognostic Model

Estimation of Immune Cell Type Fractions and estimate score

Statistical Analysis

Results

Identification of DEGs related to immune

Weighted Co-expression Network Construction And Key Modules Identification

Establishment of the lasso cox-based prognostic gene signature

Validation Of The Immuneriskscore In Tcga Crc Cohort

Correlations between the ImmuneRiskScore and immune biomarkers

Discussion

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1