Prognostic values of immune scores and immune microenvironment-related genes for hepatocellular carcinoma

Relationship of ESTIMATE scores with immune infiltration

A flowchart of the analysis procedure for this study is shown in Figure 1. We used four different methods to analyse the correlation of scores of all immune-related cell types. As shown in Figure 2A, the mean correlation of different immune cells was larger than 0.5. The ten most correlated immune cell scores with other scores were LCK (R=0.69), Co_stimulation (R=0.62), dendritic (R=0.62), Tfh (R=0.61), Co_inhibition (R=0.61), cytolytic (R=0.6), CD8_Tcell (R=0.59), ImmuneScore (R=0.59), ESTIMATEScore (R=0.58), and cytotoxic lymphocytes (R=0.57). The concentrations of immune-related cells calculated with different methods had a certain consistency. In hierarchical clustering heat maps of various scores (Figure 2B), we found immune cell scores in each sample by different methods were also consistent.

Figure 1. Flowchart describing the procedure of analyzing and validating the prognostic values of immune scores and immune microenvironment-related genes.

Figure 2. Correlations between three ESTIMATE scores and other types of immune-related scores. (A) Clustering heat map analyzed Spearman’s rank correlation coefficient. (B) Hierarchical clustering heat map using correlation to calculate distance. (C) Mean correlations of four methods using to calculating immune scores.

We further tested the average correlations between the immune scores calculated by the four methods and other types of scores (Figure 2C). It could be seen that immune scores calculated by ESTIMATE that were larger than 0.54 had the highest average correlation than that using the other three methods. This indicated that ImmuneScore, StromalScore, and ESTIMATEScore calculated by the ESTIMATE method were closely related to components of immune cells in the tumour microenvironment.

Relationship between ESTIMATE immune scores and HBV/HCV/molecular subtypes

Based on the three scores generated by the ESTIMATE algorithm, we analysed the relationship between immune scores and HBV/HCV/molecular subtypes that had been reported in previous comprehensive genomic analysis of liver cancer [13] (Supplementary file 1). As shown in Figure 3, we could see that HBV and HCV factors had no significant effect on ImmuneScore, StromalScore, and ESTIMATEScore (all P>0.05). However, the three scores in the three molecular types of expression profiles were significantly different, and iCluster3 had the lowest score (all P<0.001), suggesting a difference in the immune scores among the molecular subtypes in liver cancer.

Figure 3. Differences of the three ESTIMATE immune scores respectively in (A) HCV-positive and -negative patients, (B) HBV-positive and -negative patients, (C) three molecular subtypes based on mRNA expression profiles.

Recurrent analysis of immune scores

To observe the relationship between ESTIMATE scores and recurrence, we classified all the samples into high- and low-score subgroups according to the median ESTIMATE immune score, then the Kaplan-Meier method was used for recurrent difference analysis. As shown in Figure 4A–4C, patients in the high-score group had a higher 4-year recurrence-free rate than the low-score group, especially for progression-free 4-year survival, indicating that the immune scores were potential predictive recurrent markers.

Figure 4. Kaplan-Meier curve for recurrent analysis of overall survival time and progression-free survival time by immune scores. (A) ESTIMATEScore, (B) ImmuneScore, (C) StromalScore. H (red solid line) and L (blue solid line) respectively represented high-score group and low-score group. The red dotted line and blue dotted line respectively represented upper and lower limit of 95% confidence intervals. The vertical dotted line was at 4 years. OS, overall survival time; PFS, progression-free survival time.

Relationship between ESTIMATE immune scores and TP53, CTNNB1, and AXIN1 mutations

Many previous reports indicated that TP53, CTNNB1 (β-catenin), and AXIN1 mutations are closely related to liver cancer development. Therefore, we analysed the relationship between mutations of these three genes and ESTIMATE’s immune scores. Firstly, the mutation data of TP53, CTNNB1, and AXIN1 were extracted from SNP data treated by Mutect in TCGA. Then, the relationship between the immune scores of the mutant and non-mutant group divided by the three genes was analysed separately (Figure 5A–5C). It could be seen that StromalScore had significant differences among all three of the mutated genes, with the mutation group being smaller than wild-type group (p=0.001 for TP53, p<0.001 for CTNNB1, p=0.005 for AXIN1). ImmuneScore was significantly lower in the CTNNB1 mutant group than in the wild-type group, while ESTIMATEScores were significantly lower in the CTNNB1 and AXIN1 mutant groups than in the wild-type group.

Figure 5. (A–C) Relationship between mutations of TP53, CTNNB1, AXIN1 and ESTIMATE’s immune scores. Mut, mutation group. WT, wild type group.

Immune score-related module analysis

As shown in Supplementary Figure A, clustering analysis of samples was performed. Three hundred and ninety-six samples were finally obtained after excluding the outliers. To ensure it was a scale-free network, we chose β=12 (Supplementary Figure B, C). Finally, a total of 15 modules was obtained (Supplementary Figure D). The grey module was a gene collection that could not be aggregated into other modules. The transcripts statistics of each module are shown in Table 1. It could be seen that 6,226 transcripts were assigned to 14 co-expression modules. The correlations between 15 module eigenvectors and the three immune scores were calculated (Supplementary Figure E). The blue and yellow modules had the highest correlation with the three immune scores, and the average correlation coefficients were 0.59 and 0.48, respectively. There were 319 and 142 genes in these two modules, totalling 461 genes.

The functions of genes in the modules most relevant to the immune score were enrichment analysed (Supplementary file 2). The blue module was enriched in 85 KEGG pathways, 712 biological processes, 75 cell components, and 53 molecular functions. Meanwhile, the yellow module was enriched in 61 KEGG pathways, 280 biological processes, 35 cellular components, and 29 molecular functions. The most significant of the top 20 KEGG pathways and GO terms of the blue module are shown in Figure 6A–6D. We determined that: the KEGG pathways were mainly related to inflammation, biological processes were mainly related to T cell activation, cell compositions were mainly enriched in MHC class II protein complex, MHC protein complex, and other components, and molecular functions were mainly enriched in cytokine receptor activity, MHC class II receptor activity, and other functions. These pathways and GO terms were closely related to the immune process.

Figure 6. Enrichment analysis of genes in modules most relevant to immune scores. Top 20 (A) KEGG pathway, (B) GO BP, (C) GO CC, and (D) GO MF enriched by Blue module. Top 20 (E) KEGG pathway, (F) GO BP, (G) GO CC, and (H) GO MF enriched by Yellow module. (I) Intersections of KEGG pathway and GO Term enriched by Blue and Yellow modules. BP, biological processe. CC, cell component. MF, molecular function.

The top 20 KEGG pathways and GO terms of the yellow module are shown in Figure 6E–6H. It could be seen that: the module-enriched KEGG pathways were mainly related to focal adhesion and the PI3K–AKT signalling pathway, biological processes were mainly related to angiogenesis process, cell compositions were enriched in proteinaceous extracellular matrix, cell–substrate junctions, and other components, and molecular functions were mainly enriched in functions such as growth factor binding and transmembrane receptor protein kinase activity. These pathways and GO terms were closely related to tumorigenesis and metastasis.

Further analysis of the relationship between the two modules enriched in KEGG pathways and GO terms is shown in Figure 6I. The blue and yellow modules had only a small number of common intersections, which suggested that genes in the two modules performed different functions.

Module-related hub gene screening

To select genes related to immune score, we calculated the weight of the relationship between genes in the two modules. We selected an inter-gene weight threshold greater than 0.2, and finally obtained a weight co-expression network of the genes in the two modules. As shown in Figure 7A, the network consisted of 155 nodes and 1,475 edges, including 150 blue module genes and five yellow module genes. It could be seen that more relevant the module was, the closer the genes were to other genes in the network.

The largest network module had 36 genes, and the degree distribution of the network was further analysed. As shown in Figure 7B, the degree of most nodes was small, and a small number of nodes were large, which was consistent with the characteristics of biological networks. We calculated the correlation coefficient between genes and modules in the two co-expression modules (Figure 7C). It was more than 0.7 for most genes, indicating that genes in the module had high expression similarity. We selected 36 genes as immune microenvironment-related genes, which belonged to the largest network module, and a correlation coefficient greater than 0.7 (Supplementary file 3). Ten LCK metagenes and 26 blue module genes were included.

Figure 7. (A) Weight co-expression network of genes in Blue and Yellow modules. (B) The degree distribution of the network. (C) The correlation coefficient between genes and modules.

Network interaction analysis of immune microenvironment genes

As shown in Figure 8A, a total of 74 edges and 25 nodes were obtained. Fourteen (93.3%) of the genes interacted directly with LCK metagenes in the PPI network, except SASH3. We further compared genes that interacted directly with LCK metagenes in the PPI network with co-expressed genes (Figure 8B). Finally, we screened 12 genes that were co-expressed with immune scores and did not interact directly with LCK metagenes, which might be new genes associated with the immune microenvironment (Table 2).

Figure 8. (A) PPI network of the immune microenvironment genes. Red circle was LCK Metagenes. Light blue circle represented genes belonged to Blue module. (B) Venn diagram showing the number of genes in WGCNA, StringPPI and LCK Metagenes. WGCNA was from co-expression network. StringPPI was from direct interaction with LCK Metagenes in protein interaction network. MetaGenes was a merger of 13 set of MetaGenes.

Relationship between 12 new immune microenvironment genes and immune checkpoint genes

We further analysed the correlation between the 12 new immune microenvironment genes and eight immune checkpoint genes, including

PDCD1, CD274, PDCD1LG2, CTLA4, CD86, CD80, CD276, and VTCN1. As shown in Figure 9, there were significantly positive correlations between the new immune microenvironment genes and seven immune checkpoint genes, except VTCN1. The average correlation coefficient was greater than 0.52, suggesting that these immune microenvironment genes might be potential targets for immunotherapy.

Figure 9. Correlation between the 12 immune microenvironment genes (blue green bar) and 8 immune checkpoint genes (light green bar). Lower panel was Pearson correlation coefficient between each gene. Upper panel was scatter plot of expressions between each gene. Diag panel was expression of each gene.

Prognostic value of 12 immune microenvironment genes

From Figure 10, we could see that, except for RCSD1, high expression of the other 11 genes was significantly associated with a good prognosis (P<0.05).

Figure 10. Kaplan-Meier curve for prognostic analysis of the 12 immune microenvironment genes. The red dotted line and blue dotted line respectively represented upper and lower limit of 95% confidence intervals of gene expression.

External validation of the expression of the 12 immune-related genes and their relationship with immune scores

In the GSE62232 dataset, two immune microenvironment gene expression profiles were not detected, so the other 10 genes were used to calculate their correlation with immune scores. As shown in Figure 11A, nine genes, except TMC8, had significant correlations with the immune score (P<0.001), which was consistent with the analysis in TCGA cohort.

Figure 11. External validation of 12 immune-related genes. (A) Pearson correlation of expression and ImmuneScore using GSE62232 dataset. TRBC2 and TRBV28 were not detected. 9 genes, except TMC8, had significant correlation. Y-axis represented expression level of the genes in each sample. (B) Representative immunohistochemical staining images of 7 genes in normal liver tissue and HCC specimen. Images were taken from the Human Protein Atlas (https://www.proteinatlas.org). TMC8 was strongly positive, BIN2, GIMAP7, SPOCK2 were moderately positive, and FYB1, RCSD1, SASH3 were weakly positive in HCC tissues relative to their expression levels in normal liver tissues.

In the validation of protein expression encoded by the genes using clinical specimens from Human Protein Profiles, we found that TMC8 was strongly positive, BIN2, GIMAP7, and SPOCK2 were moderately positive, and FYB1, RCSD1, and SASH3 were weakly positive in HCC tissues relative to their expression levels in normal liver tissues (Figure 11B). However, NAPSB, TRAC, TRBC2, GPSM3, and TRBV28 were not found on the website.

Validation of expression of 3 immune-related genes in HCC tissues

To further validate the results, qRT-PCR was applied to analyse the relative mRNA expression of TMC8 and BIN2 in fresh HCC and adjacent non-tumour tissues. The results showed that elevated expression of TMC8 and BIN2 relative mRNA in tumor tissues compared to adjacent normal tissues (Figure 12A–12B). Our western blotting analysis demonstrated that HCC tissues exhibited relative higher levels of TMC8,BIN2 and SPOCK2 protein expression than those in normal liver tissues (Figure 12C–12F). In addition, IHC analysis was also conducted to determine TMC8,BIN2 and SPOCK2 protein expression level in HCC. From the immunostaining, we could observe the positive immunoreactivity were primarily localized in the cytoplasm (Figure 12G). The three proteins were all up-regulated in tumor tissues compared with the adjacent normal liver tissues (Figure 12H–12J). These results all indicated that TMC8,BIN2 and SPOCK2 were overexpressed in HCC tissues.

Figure 12. Measurement of TMC8, BIN2 and SPOCK2 at mRNA and protein level in our cohort. (A, B) Relative mRNA levels of TMC8 and BIN2 in 10 HCC samples were both overexpressed compared with matched normal samples by qRT-PCR. (C) Representative western blotting images showed protein expression of TMC8, BIN2 and SPOCK2 were overexpressed in HCC tissue than those in normal liver tissue. (D–F) Western blotting analysis demonstrated that mean greyscale of TMC8, BIN2 and SPOCK2 were all higher in 10 fresh-frozen HCC tissues compared with those in matched adjacent normal liver tissues. (G) Representative immunohistochemical staining images of TMC8, BIN2 and SPOCK2, which were all mainly expressed in the cytoplasm. Original magnification: x200. (H–J) Mean protein expression of TMC8, BIN2 and SPOCK2 in 35 HCC were all significantly higher compared with those in adjacent non-tumour tissue by immunohistochemistry. MAD: mean areal density.

Data from TCGA cohort

The transcriptome expression profiles and corresponding clinical information of HCC were downloaded from the Genomic Data Commons Application Programming Interface of TCGA (https://cancergenome.nih.gov/). Surgical resection samples were collected from patients who were diagnosed with HCC and did not receive prior adjuvant treatment for their disease. FPKM (fragments per kilobase million) data of RNA-Seq was first downloaded and then transferred to TPM (transcripts per million) expression data. The expression data also included single nucleotide polymorphism (SNP) expression data, containing 371 cases and 424 files as of February 2019.

Data analysis

Four different algorithms were used to calculated the scores of related immune cells. The expression levels of 13 immune metagenes were listed, including ImmuneScore, LCK (lymphocyte-specific protein tyrosine kinase), TFH, Tregs, cytolytic, MHC1, MHC2, NK, macrophages, STAT1, IF_I, Co_stimulation, and Co_inhibition, which corresponded to various immune cell types and reflected various immune functions [10] (Supplementary file 2). Based on the gene expression levels of these metagenes, their scores were calculated using the median expression level of each gene in all the samples (Supplementary file 3).

There were large numbers of immune cell-specific genes highly expressed in the tumor microenvironment. The scores of six immune cell types of liver cancer, including B cells, CD4 T cells, CD8 T cells, neutrophils, macrophages, and dendritic cells, were downloaded from the Timer (https://cistrome.shinyapps.io/timer/) database, which used constrained least squares fitting on the informative immune signature genes to detect the immune cell infiltration in tumour tissue [27] (Supplementary file 4).

ESTIMATE is an algorithmic tool for predicting tumour purity, which uses the gene expression profiles of 141 immune genes and 141 stromal genes to generate ESTIMATE scores [9]. The presence of infiltrated immune cells and stromal cells in tumour tissues were calculated using related gene expression matrix data, represented by ImmuneScore and StromalScore, respectively (Supplementary file 5).

The abundance of eight immune cell types, including T cells, CD8 T cells, cytotoxic lymphocytes, NK cells, B lineage, monocytic lineage, myeloid dendritic cells, neutrophils, and other two types of cells, including endothelial cells and fibroblasts, in each sample were estimated with MCPcounter software (Supplementary file 6).

To explore the relationship of ESTIMATE scores with other immune scores, we analysed the scores of several immune-related cell types using the four different algorithms.

Relationship between ESTIMATE’s immune score and TP53, CTNNB1, and AXIN1 mutations was also analysed.

Mining immune score-related modules

To further explore prognostic markers related to the immune microenvironment of liver cancer, we obtained 423 expression profiles of all samples. More than 75% transcripts, totalling 11,803 TPM >1 and median absolute deviation greater than the median were selected in these samples. Firstly, hierarchical clustering analysis of the samples was undertaken, and outliers with a height over 100,000 were excluded. Distances between each transcript were calculated with Pearson’s correlation coefficient. WGCNA was constructed using the R package WGCNA.

The co-expression network conformed to a scale-free network, which was to say log(k) of the node with connection degree k was negatively correlated with log(P(k)) of probability of the node, and the correlation coefficient was greater than 0.85. To ensure it was a scale-free network, we chose β=12 (Supplementary Figure). Then, the expression matrix was converted into an adjacency matrix, which was converted into a topological matrix. Hierarchical clustering analysis of genes was performed by the average linkage hierarchical clustering method. According to criteria of the hybrid dynamicTreecut software package, each lncRNA network module had a minimum of 30 genes. After confirming the gene modules, we calculated the eigengenes of each module in turn, then clustered the modules, and merged the closer modules into new ones by setting height=0.25, deepSplit=2, and minModuleSize=30. We further analysed the functions of genes that were most correlated with the immune score in the modules. KEGG and GO enrichment analysis of the genes were performed with FDR <0.05 using R package clusterProfiler.

Constructing the PPI network

We analysed PPI network of these 36 immune microenvironment-related genes using R package. We mapped the 36 genes from STRING database, selected interaction scores greater than 0, and then obtained the network relationship among these genes.

Prognostic value of 12 immune microenvironment-related genes

To analyse the prognostic value of 12 immune microenvironment-related genes, all samples were divided into high and low expression groups according to their expression profiles. Then, prognostic differences between the two groups were analysed by the K-M method.

External validation in GEO datasets and Human Protein Profiles

To ensure the accuracy of results from TCGA cohort, Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) datasets were analysed to verify the relationship between the 12 genes and immune scores. GSE62232, containing 81 liver cancer samples and 10 normal samples, was selected. The standardized expression matrix and sample information files were downloaded, and immune scores for each sample were calculated using R ESTIMATE package. We also extracted the gene expression profiles of the 12 immune-related genes, and then calculated Pearson’s correlation coefficient between the genes and the immune scores.

To further verify the accuracy of the results, we analysed the expression of the proteins encoded by the 12 immune-related genes using clinical specimens from The Human Protein Atlas (https://www.proteinatlas.org) [28].

The data from TCGA, GEO and The Human Protein Atlas were all publicly available and open access, so no approval was needed from the ethics committees.

Quantitative real-time polymerase chain reaction (qRT-PCR) evaluation

From January 2017 to December 2018, 10 HCC patients who underwent curative resection and didn’t receive neoadjuvant therapy before surgery at the First Affiliated Hospital of Zhengzhou University (Zhengzhou, China) participated in this study in accordance with the provisions of Helsinki Declaration.

To evaluate the gene expression of TMC8 and BIN2, 10 pairs of fresh HCC and adjacent normal samples were collected and tested using qRT-PCR method. Total RNA was extracted with Trizol reagent (Invitrogen, USA) and reverse-transcribed into cDNA using RevertAid First Strand cDNA Synthesis Kit (Thermo, USA) according to the manufacturer’s instruction. qRT-PCR was performed with the FastStart Universal SYBR Green Master (Rox) (Servicebio, China). The quantitative levels of gene expressions in HCC and adjacent paracancerous tissue were evaluated using the 2-ΔΔCt relative quantification method. The primers were as follows: TMC8: forward primer, 5′-GACTCTGCTGGGTCAGGGCTAT-3′, reverse primer, 5′-TCCACCTTGAACTCGTTGCTG-3′; BIN2: forward primer, 5′-ACGAGGAGAAACTGGCTGACC-3′, reverse primer, 5′-CACTGTCATAGTCCACGAGTTTCC-3′; β-actin (used as control): forward primer, 5′-GGAAGCTTGTCATCAATGGAAATC-3′, reverse primer, 5′-TGATGACCCTTTTGGCTCCC-3′.

Western blotting analysis

Proteins were isolated from these 10 pairs of fresh-frozen HCC and adjacent normal tissues as previously described [29]. Protein concentrations of TMC8, BIN2 and SPOCK2 were measured using bicinchoninic acid assay. Equal amounts of proteins were separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to a polyvinylidene difluoride (PVDF) membrane (Millipore, USA). After blocking with 5% skimmed milk, the membrane was incubated with primary rabbit polyclonal antibodies (TMC8, BS-13116R, 1:500, Bioss, China; BIN2, BS-9727R, 1:500, Bioss, China; SPOCK2, AB217044, 1:1000, Abcam, UK) and ACTIN (GB11001, 1:3000, Servicebio, China) overnight at 4 °C, respectively. After washing the membrane, it was then incubated with a horseradish peroxidase-conjugated goat anti-rabbit secondary antibody (GB23303, 1:3000, Servicebio, China). Finally, the bound immunocomplexes were detected using ECL methods.

Immunohistochemical (IHC) validation of the immune-related genes

All haematoxylin and eosin stained slides of tumor and adjacent non-tumour samples from another 35 HCC patients, confirmed by two experienced pathologists, were collected. Three immune-related genes, TMC8, BIN2 and SPOCK2, were selected to be further validated with immunohistochemical method. All formalin-fixed paraffin-embedded HCC samples were collected to examine the protein levels of the three immune-related genes.

Slides with 4um sections from the paraffin-embedded specimens were deparaffinized and rehydrated. After conducting heat-mediated antigen retrieval, endogenous peroxidase activity was blocked by incubating the sections with 3% hydrogen peroxide at room temperature for 25 min. Then the slides were washed with phosphate buffered saline (PBS) three times and incubated with bovine serum. They were further incubated with primary and secondary antibodies. After washing with PBS and incubating with 3.3′-diaminobenzidine solution for 5 min, each section was counterstained with hematoxylin. All IHC slides were analyzed using Image-pro plus 6.0 software (Media Cybernetics, USA). The results were tested by two pathologists. For TMC8, BIN2 and SPOCK2 expression analysis, the relevant primary antibodies (BS-13116R, 1:400, Bioss, China; BS-9727R, 1:200, Bioss, China and AB217044, 1:200, Abcam, UK) were used. The results were all expressed with mean areal density (MAD).

Statistical methods

All analyses were conducted using R software (version 3.5.2). ImmuneScore, StromalScore and ESTIMATEScore of each sample were calculated with R package ESTIMATE. The abundance of eight immune cell types, endothelial cells, and fibroblasts in each sample were estimated with MCPcounter software. Distances between each transcript were calculated with Pearson’s correlation coefficient. WGCNA was constructed using the R package WGCNA. The relationships between immune scores and HBV/HCV/molecular subtypes were tested by Analysis of Variance test. Survival R package was used for Kaplan-Meier curve analysis. Wilcox test was performed to analyze relationships between ESTIMATE immune scores and TP53, CTNNB1, AXIN1 mutations. Correlations between immune-related genes and immune checkpoint gene, ImmuneScores were studied by corrgram R package. The two-tailed paired t-test was used for data of TMC8, BIN2 and SPOCK2 expression at mRNA and protein level. All statistical tests were two sided and p values <0.05 were considered as statistically significant.