Identification of Key Functional Modules and Immunomodulatory Regulators of Hepatocellular Carcinoma.

Despite the advances in the treatment of hepatocellular carcinoma (HCC), the prognosis of HCC patients remains unsatisfactory due to postsurgical recurrence and treatment resistance. Therefore, it is important to reveal the mechanisms underlying HCC and identify potential therapeutic targets against HCC, which could facilitate the development of novel therapies. Based on 12 HCC samples and 12 paired paracancerous normal tissues, we identified differentially expressed mRNAs and lncRNAs using the "limma" package in R software. Moreover, we used the weighted gene coexpression network analysis (WGCNA) to analyze the expression data and screened hub genes. Furthermore, we performed pathway enrichment analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. In addition, the relative abundance of a given gene set was estimated by single-sample Gene Set Enrichment Analysis. We identified 687 differentially expressed mRNAs and 260 differentially expressed lncRNAs. A total of 6 modules were revealed by WGCNA, and MT1M and MT1E genes from the red module were identified as hub genes. Moreover, pathway analysis revealed the top 10 enriched KEGG pathways of upregulated or downregulated genes. Additionally, we also found that CD58 might act as an immune checkpoint gene in HCC via PD1/CTLA4 pathways and regulate the levels of tumor-infiltrating immune cells in HCC tissues, which might be an immunotherapeutic target in HCC. Our research identified key functional modules and immunomodulatory regulators for HCC, which might offer novel diagnostic biomarkers and/or therapeutic targets for cancer immunotherapy.

1. Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignant tumors and the fourth leading cause of cancer-related death globally, which has posed a substantial financial and health burden [1]. In China, HCC accounts for about 90% of all cases of primary liver cancer, and 422,100 individuals succumb to liver cancer annually [2]. Despite the advances in HCC treatment (especially immunotherapy), its mortality holds an increasing trend and the prognosis of patients with HCC remains unsatisfactory. However, the complex mechanisms underlying the initiation and progression of HCC make it challenging to develop novel therapeutic strategies.

Accumulating studies reveal the critical roles of immune cells in the initiation, metastasis, and recurrence of HCC [3, 4]. Multiple regulatory molecules could inhibit the antitumor activity of tumor-associated immune cells, thus resulting in immune escape [5-7]. As a strategy to normalize the antitumor immune responses against cancer cells, cancer immunotherapy has achieved clinical success in the last five years and revolutionized the treatment landscape of HCC [8]. However, the limited response and significant toxicity impede the therapeutic effect. Therefore, it is necessary to reveal the underlying mechanisms of HCC and identify potential therapeutic targets, with an emphasis on immunomodulatory regulators.

This study identified differentially expressed genes and lncRNAs based on the expression profile of 12 HCC samples and 12 paired paracancerous normal tissues. The weighted gene coexpression network analysis (WGCNA) was used to screen vital modules and hub genes correlated to HCC. Moreover, we performed pathway enrichment analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. In addition, we estimated the relative abundance of a given gene set by single-sample Gene Set Enrichment Analysis (ssGSEA). This study is aimed at identifying key functional modules and immunomodulatory regulators for HCC, which might offer novel diagnostic biomarkers and/or therapeutic targets for cancer immunotherapy.

2. Materials and Methods

2.1. Data Acquisition and Preprocessing

The expression profile of GSE115018 was downloaded from the Gene Expression Omnibus database. In GSE115018, Shi et al. [9] collected 12 HCC samples and 12 paired paracancerous normal tissues from the first affiliated hospital of Guangxi Medical University. All the patients received the first operation on primary disease without radiochemotherapy. The diagnosis of HCC was further confirmed using histopathology after surgery. Within 30 minutes after isolation, the tissues were immediately frozen in liquid nitrogen and stored at -80°C. The Ethics Committee of the First Affiliated Hospital of Guangxi Medical University approved this study, and informed consent for participation was acquired from all participants and their families [9]. We matched probes with gene symbols after removing redundant data (time, null value, etc.). A total of 9,945 mRNAs and 5,072 lncRNAs were analyzed.

2.2. Differential Expression Analysis

We first used the “wateRmelon” package to correct background, normalize quantile, and summarize quantile to eliminate potential error. Then, we used the “limma” package to analyze the differentially expressed mRNAs and lncRNAs between HCC and paired paracancerous normal tissues with criteria of adjusted P value < 0.05 and ∣log2 fold − change (FC) | ≥1.

2.3. WGCNA and Screening of Hub Genes

Weighted gene coexpression network analysis (WGCNA) is a systematic biology method [10], which identifies highly correlated genes and related modules to external sample traits. Before network construction, we removed obvious outlier samples or samples with excessive numbers of missing entries to avoid data deviation. Then, the step-by-step network construction and module detection were performed [10].

According to the criterion of approximate scale-free topology, we chose the soft-threshold power (β-value) of 7 to determine a scale-free topology index (R2) of 0.86. Then, we calculated adjacencies using the soft-threshold power and transformed the adjacency into the Topological Overlap Matrix (TOM), which was an average linkage hierarchical clustering with a dissimilarity measure to detect gene modules. The minimum module size of 30 and cut height of 0.99 were set to identify and merge gene modules with similar expression probes.

Furthermore, we used Zsummary to evaluate the functional module preservation [11]. Zsummary is comprised of four statistics related to density and three connectivity-related statistics, and it was created to quantitatively assess whether the density and connectivity patterns of modules defined in a reference dataset are preserved in a test dataset. A Zsummary value between 2 and 10 indicates moderate module preservation, whereas a Zsummary > 10 provides strong support for module preservation.

Genes with the highest degree of connectivity in a module were identified as hub genes, which could determine the biological significance of the module. We correlated the different module eigengenes (MEs) and the clinical traits. The gene significance (GS) quantifies the association of individual genes with the clinically interesting trait, and the module membership (MM) acts as the correlation between MEs and the gene expression profiles. Hub genes would be chosen if GS value was >-log10 (0.05) and the absolute value of kME was >0.85.

2.4. Pathway Enrichment Analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological system from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies [12]. The logFC of DEGs obtained from differential expression analysis was applied for enrichment analysis. R packages of “clusterProfiler” and “enrichplot” were used to perform KEGG enrichment analysis with a threshold of P value <0.05. Further, enriched signaling pathways were visualized with “dotplot” and “gseaplot” packages of R software.

2.5. Network Analysis and Target Relationship Prediction

Coexpression network analysis was conducted based on the Pearson correlation coefficient to explore the correlations between differentially expressed lncRNAs and mRNAs. R was applied to the output node and edge files of the eligible paired lncRNA-mRNA interaction, and Cytoscape 3.6.0 software was used to visualize the lncRNA-mRNA coexpression network.

Then, miRNA targets of lncRNAs and miRNAs were carried out based on the miRanda database (http://www.microrna.org/) and TargetScan 7.2 (http://www.targetscan.org/vert_72/). The miRNAs cotargeting the certain lncRNA and the certain mRNA were filtered using Venny 2.1.0 (http://bioinfogp.cnb.csic.es/tools/venny/).

2.6. Expression Analysis of Hub Genes Based on TCGA

The expression of hub genes in HCC and normal tissues was evaluated using the UALCAN, an interactive web portal to perform in-depth analyses of TCGA gene expression data [13]. Meanwhile, the expression levels of hub genes were also evaluated in patients with various tumor grades, ages, genders, and races.

2.7. TIMER Database Analysis

We used the TIMER database to analyze the association between the hub gene expression and the abundance of infiltrating immune cells, including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells [14].

2.8. The Immune Cell Infiltration in Tumor Cells

We comprehensively estimated the infiltration level of immune cell populations of each sample by ssGSEA, which was implemented in the GSVA package [15]. As an extended gene set enrichment analysis method, ssGSEA was designed to compute the separate enrichment scores for a particular gene set in each sample instead of the gene-phenotype association score. The gene set from a previous study [16], which included 364 genes representing 24 microenvironment cell types (see their Supplementary Table S1), was input into the ssGSEA algorithm. The overall immune cell infiltrating levels were estimated by the R estimate package [17].

3. Result

3.1. DEG Identification

After correcting background, normalizing quantile, and summarizing quantile (Figures 1(a) and 1(b)), we identified a total of 687 differentially expressed mRNAs in HCC tissues compared with paracancerous normal tissues. Figure 1(c) shows the volcano plot of differentially expressed mRNAs. The top 10 differentially upregulated genes (including SPINK1, HIST1H2AG, HIST1H3H, CCNA2, KIF20A, MKI67, MDK, CCNE2, GPC3, and TRIM31) and top 10 downregulated genes (including CYP2E1, MT1M, MT1E, SLC25A47, MT1G, COLEC10, FCN3, CRHBP, DCN, and SAA1) are shown in the heat map (Figure 1(d)).

Figure 1

Identification of differentially expressed mRNAs. Box plots show the distribution of the relative mRNA expression in each sample before (a) and after (b) normalization of GSE115018. Each box corresponds to one sample. The middle line corresponds to the median. (c) Volcano plot of differentially expressed mRNAs. Differentially expressed mRNAs were screened with criteria of adjusted P value < 0.05 and ∣log2 fold − change (FC) | ≥1. (d) The cluster heat map of the top 10 upregulated and downregulated mRNAs.

Moreover, a total of 260 differentially expressed lncRNAs were identified from dataset GSE115018 after preprocessing (Figures 2(a) and 2(b)). As shown in the volcano plot (Figure 2(c)), 102 of these lncRNAs were upregulated in the tumor, while 158 of them were downregulated. The top 10 upregulated and downregulated lncRNAs in HCC tissues compared with paracancerous normal tissues are shown in the heat map (Figure 2(d)).

Figure 2

Identification of differentially expressed lncRNAs. Box plots show the distribution of the relative lncRNA expression in each sample before (a) and after (b) normalization of GSE115018. (c) Volcano plot of differentially expressed lncRNAs. Differentially expressed lncRNAs were screened with criteria of adjusted P value < 0.05 and ∣log2 fold − change (FC) | ≥1. (d) The cluster heat map of the top 10 upregulated and downregulated lncRNAs.

3.2. WGCNA and Screening of Hub Genes

To identify key modules related to HCC, the soft-threshold power of 7 was set to ensure a scale-free network (scale R2 = 0.86) (Figure 3(a)). With a minModuleSize of 30 and a CutHeight of 0.99, a total of six modules were identified in WGCNA: blue, brown, green, red, turquoise, and yellow module (Figure 3(b)). The heat map depicting the TOM of genes is shown in Figure 3(c). The darker parts indicated a higher degree of connectivity and suggested that these genes might be highly related to HCC. The results of GS showed that the module significance (MS) of the turquoise and red modules was higher than that of the other modules. The module eigengene tree showed a similarity among different modules (Figure 3(d)). The module yellow and module blue had high similarity, and the correlation coefficient between these two modules was 0.75. Scatterplots of module membership and GS showed where each gene was in the relevant network (Figures 3(e) and 3(f)). Additionally, we also performed enrichment analysis on the blue and yellow modules using Gene Ontology (GO) and the KEGG database. As shown in SFigure 1A, the enriched biological processes were mainly involved in fatty acid metabolic process, response to xenobiotic stimulus, cellular response to xenobiotic stimulus, xenobiotic metabolic process, and drug catabolic process. The cellular components were primarily enriched in mitochondrial matrix and mitochondrial inner membrane, whereas the enriched molecular functions were mainly associated with steroid hydroxylase activity and monooxygenase activity. KEGG pathway analysis suggested that the metabolism of xenobiotics by cytochrome P450 was the most enriched pathway, followed by complement and coagulation cascades, glucagon signaling pathway, lipid and atherosclerosis, and PI3K-Akt signaling pathway (SFigure 1B).

Figure 3

Construction of weighted coexpression network and module analysis. (a) Soft-threshold selection. A soft-threshold power of 7 with an R2 of 0.87 was selected. (b) Cluster dendrogram. Each color indicates one specific coexpression module, and each vertical line indicated an individual gene. Six modules were identified in WGCNA, including blue, brown, green, red, turquoise, and yellow module. (c) Network heat map plot of all genes suggesting the interaction of coexpression genes via TOM dissimilarity. The axe colors indicate the respective modules. The color intensity represents the degree of overlap, where a darker yellow suggests a higher degree of connectivity. (d) Barplot of mean gene significance associated with HCC. A higher value shows a more significant relationship. (e) Relation between different module eigengenes. (f) Scatterplots of genes in the six modules. Module membership is set as the x-axis, whereas the y-axis is gene significance.

Zsummary of brown, green, red, turquoise, and yellow modules were all <10 (Figure 4(a)), and we chose these five modules to perform further analysis and identify hub genes. With a threshold of GS > −log10 (0.05) and the absolute value of kME > 0.85, we selected hub genes and intersected with the top 10 upregulated or downregulated ones. Hub genes in the five modules are shown in Figure 4(b). Importantly, MT1M and MT1E genes from the red module were identified as hub genes. Moreover, we constructed a network for genes from the red module, and MT1E and MT1M genes were in the network center (Figure 4(c)).

Figure 4

The Zsummary statistics of the module preservation and hub gene screening. (a) Zsummary score analysis of different modules. The dashed blue and green lines suggested the thresholds of 2 and 10, respectively. A Zsummary value between 2 and 10 indicates moderate module preservation, whereas a Zsummary > 10 provides strong support for module preservation. (b) Hub genes in different modules. The red solid points present the hub genes in each module. (c) Network of genes in the red module.

3.3. KEGG Pathway Analysis

We performed KEGG pathway analysis to achieve a more in-depth insight into the biological roles of the identified DEGs. Figures 5(a) and 5(b) show the top 10 enriched KEGG pathways of upregulated and downregulated genes. For the upregulated genes, complement and coagulation cascades, retinol metabolism, and glycolysis were the top three pathways. Moreover, DNA replication, alcoholism, and viral carcinogenesis were the top three enriched pathways for downregulated genes. Different genes enriched by these pathways were shown in Figure 5(c). Notably, the retinol metabolism pathway was highly expressed in HCC, and the coexpression network of lncRNAs and mRNAs from the retinol metabolism pathway was shown in Figure 5(d).

Figure 5

KEGG pathway analysis and coexpression network. (a) Dot plot of the top 10 enriched KEGG pathways of upregulated or downregulated genes. The logFC of DEGs obtained from differential expression analysis was applied for enrichment analysis. (b) The mountain range of the top 10 enriched KEGG pathways of upregulated or downregulated genes. (c) Cnetplot of the correlation between metabolism pathways and differentially expressed mRNAs. (d) Coexpression of lncRNA and mRNA network from the retinol metabolism pathway.

3.4. Expression Analysis of MT1E and MT1M Based on TCGA

Figures 6(a) and 6(d) show an overview on the MT1E and MT1M expression levels in tumors compared with normal samples across multiple cancer types in the TCGA database. MT1E and MT1M were significantly downregulated in HCC tissues than in normal samples (both P value < 0.01; Figures 6(b) and 6(e)). Moreover, Figures 6(c) and 6(f) showed that this trend was consistent regardless of tumor grade (grades 1 to 4), age (21-40, 41-60, 61-80, and 81-100 years old), gender (male and female), and race (Caucasian, African-American, or Asian ethnicity).

Figure 6

MT1E and MT1M expression levels in HCC based on UALCAN portal analysis and their association with immune cell infiltration. (a) Boxplot showing the overall comparison of MT1E expression between tumor and normal samples across multiple cancer types in the TCGA database. (b) The expression of MT1E in normal tissues compared with HCC tissues. (c) The comparison of MT1E expression in normal tissues compared with HCC tissues from different tumor grade (grades 1 to 4), age (21-40, 41-60, 61-80, and 81-100 years old), gender (male and female), and race (Caucasian, African-American, or Asian ethnicity). (d) Boxplot showing the overall comparison of MT1M expression between tumor and normal samples across multiple cancer types in the TCGA database. (e) The expression of MT1M in normal tissues compared with HCC tissues. (f) The comparison of MT1M expression in normal tissues compared with HCC tissues from different tumor grade (grades 1 to 4), age (21-40, 41-60, 61-80, and 81-100 years old), gender (male and female), and race (Caucasian, African-American, or Asian ethnicity). (g) Correlation of MT1E and MT1M expression with immune infiltration level in HCC via the TIMER database. Both MT1E and MT1M expression were significantly correlated with the infiltration levels of B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells in HCC.

3.5. TIMER Database Analysis

We used the TIMER database to explore the correlation between the expression level of MT1E and MT1M and immune cell infiltration. Both MT1E and MT1M expression were significantly correlated with the infiltration levels of B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells in HCC. Furthermore, there was a negative correlation between MT1E expression and the infiltration of B cells, CD4+ T cells, and macrophages, while the expression of MT1M was positively associated with the infiltration of B cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells. These results indicated that MT1E and MT1M might play an important role in regulating immune infiltration in HCC (Figure 6(g).

3.6. The Association between Functional Modules and the Immunomodulatory Regulators in HCC

To further explore the biological roles of the functional modules, we conducted a correlation analysis between the functional modules and immunomodulatory regulators in the TCGA liver cancer cohort, which had a larger sample size. The relative expression abundances of the functional modules were estimated by ssGSEA (Figure 7(a)). Specifically, we observed that the green module was positively correlated with TNFSF13B, CD86, CD226, TNFSF8, CD274, and PDCD1LG2. Particularly, CD86, CD274, and PDCD1LG2 were well-recognized immune checkpoint genes involved in PD1/CTLA4 pathways. In contrast, the blue module was negatively correlated with TNFRSF18 (correlation coefficient = −0.47), TNFSF15 (correlation coefficient = −0.49), or CD58 (correlation coefficient = −0.4), whereas yellow was negatively correlated with TNFSF15 with a correlation coefficient of -0.43 (all P value < 0.05).

Figure 7

The association between functional modules and the immunomodulatory regulators in HCC. (a) The relative expression abundances of the functional modules estimated by ssGSEA. (b) The CD58 RNA expression level in tumor tissues. CD58 was highly expressed in HCC tissues with high infiltrating levels of immune cells (P value < 0.001) (c) The disease-free survival and (d) the overall survival for patients with low (below median CD58 expression value, n = 91) or high (above median CD58 expression value, n = 91) expression level. Kaplan-Meier survival curves were generated based on TCGA data using the GEPIA2 tool.

Among these immunomodulatory regulators, CD58 was highly expressed in HCC tissues with high infiltrating levels of immune cells (Figure 7(b), Wilcoxon-rank sum test, P < 0.001). Furthermore, high CD58 expression was associated with shorter disease-free survival (P value = 0.0086) and overall survival (P value = 0.0078), suggesting that CD58 was closely associated with HCC prognosis (Figures 7(c) and 7(d)).

4. Discussion

This study identified 687 differentially expressed mRNAs and 260 differentially expressed lncRNAs based on 12 HCC samples and 12 paired paracancerous normal tissues. WGCNA revealed a total of 6 modules associated with HCC, and MT1M and MT1E from the red module were identified as hub genes. Then, the pathway enrichment analysis revealed the top 10 enriched KEGG pathways, and we created a coexpression network of lncRNAs and mRNAs from the retinol metabolism pathway. Interestingly, our research indicated that CD58 might act as an immune checkpoint gene in HCC and regulate the levels of tumor-infiltrating immune cells in HCC tissues, which might be an immunotherapeutic target in HCC.

Despite the remarkable advances in the treatment strategies [18], HCC is still the fourth leading cause of cancer-related death globally, and it has posed a substantial health burden [19]. Bioinformatics technology provides a novel method to reveal the mechanisms underlying HCC and identify potential therapeutic targets, which contributes to the development of novel agents against HCC.

Metallothionein (MT) family members are reported to be associated with the cellular metabolism of metal ions and cancer development. For example, MT family members were reported to take part in the malignant transformation of hepatocytes [20]. They could also protect cells from anticancer agents and irradiation-induced damage [21]. MT1 is one of the MT family members, which is expressed in all eukaryotes. It was reported that MT1 was downregulated in HCC, and the silence of MT1 could promote the proliferation of liver cancer [22]. Moreover, MT1M was demonstrated as a tumor suppressor gene downregulated in HCC, which would contribute to liver tumorigenesis by increasing cellular NF-κB activity [23, 24]. Consistent with previous researches, our analysis discovered that MT1M and MT1E were downregulated in HCC tissue and identified them as hub genes in HCC. However, the role of MT1E in liver cancer is seldom studied. Our analysis disclosed that MT1E was downregulated in liver cancer and may serve a similar function as MT1M.

KEGG pathway analysis suggested that multiple pathways (including complement and coagulation cascades, retinol metabolism, glycolysis pathways, and so forth) were significantly enriched. Among them, the complement system and coagulation cascade are classic immunomodulatory regulators for both innate and adaptive immune responses [25]. Moreover, retinol (vitamin A) lies at metabolic crossroads of multiple biochemical reactions, which serve as a necessary precursor for retinoid and retinoic acid [26]. Immunomodulatory effects of retinoids have been revealed in various human cellular lineages, including thymocyte, lung fibroblast, Langerhans' cell, natural killer cell, peripheral blood mononuclear cell, and tumoral cell [27]. For example, retinoic acid significantly improves the expansion of Foxp3+ inducible regulatory T cell and inhibits the differentiation of T helper 17 cell, which contributes to maintaining immune homeostasis [28, 29]. In contrast, retinol deficiency increases proinflammatory cytokines and elevates T helper type 1 response [30]. The immunomodulatory roles of retinoids have been summarized in several reviews [31, 32]. The close connection between retinoids and liver diseases has also long been recognized [33, 34], and several experimental studies indicated the beneficial effects of retinoids on blocking HCC development [35-37]. Additionally, accumulating evidence suggests the immunomodulatory role of butyrate, which could affect the epigenetic status of immune cells via downregulating the enzymatic function of histone deacetylase [38-40].

In further network analysis on the retinoid metabolism pathway, CYP1A2, CYP2A13, CYP2A7, CYP2B6, CYP2C19, CYP4A11, ADH1B, ADH1C, ADH4, ADH6, and RDH5 were in the center of the network. High expression of CYP1A2 was reported to serve as a biomarker to predict recurrence-free survival of HCC [41]. Moreover, the CYP2A13, CYP2A7, CYP2B6, and CYP2C19 genes belong to the CYP2 family [42]. CYP2A13 is a human cytochrome P450 (P450) enzyme, which is widely expressed in the liver [43, 44]. Moreover, it is responsible for the metabolism of nicotine, coumarin, and tobacco-specific nitrosamine [45]. Recently, the interaction between CYP2A13 and ABCB1was reported to be closely associated with lung cancer, and CYP2A13 was identified as a potential critical metabolic enzyme gene in the carcinogenesis of lung cancer [46, 47]. Additionally, the CYP2A7 pseudogene transcript was demonstrated to affect the expression of CYP2A6 in the liver as a decoy for miR-126, but the role of CYP2A7 in HCC remains vague [48]. CYP2B6 and CYP2C19 were previously reported as unfavorable prognosis markers in breast cancer [49, 50]. In the analysis of circulating DNA from patients with advanced hepatocellular carcinoma, CYP2B6, BAX, and HNF1A genes showed the highest mutation frequency and a significant association with the clinicopathological characteristics of HCC, which suggested potential roles as driver genes in a specific clinical setting [51]. CYP2C19 belongs to cytochrome P2C subfamily members, which are known to be involved in clinical drug metabolism. The expression levels of CYP2C8, CYP2C9, and CYP2C19 genes were identified as potential prognostic markers of HCC following hepatectomy [42]. Also, the downregulation of the CYP2C19 gene is associated with aggressive tumor potential and poorer recurrence-free survival of HCC [52, 53]. In our analysis, CYP2B6 and CYP2C19 were discovered to be upregulated in HCC, and the roles in the process of liver cancer are required to be further verified.

ADH4 is an important member of the ADH family, which was involved in metabolizing a large variety of substrates such as ethanol and retinol. ADH4 mRNA and protein expression levels were markedly reduced in HCC tissues and were reported to be recognized as potential prognostic biomarkers for HCC patients. HCC patients with lower ADH4 showed a worse overall survival rate compared with those with high expression (P < 0.001), and the expression of ADH4 was an independent predictor of overall survival (HR, 0.154; 95% CI, 0.044-0.543; P = 0.004).

To further explore the biological roles of the functional modules, we also conducted a correlation analysis between the functional modules and immunomodulatory regulators. We found that green, blue, and yellow modules were closely associated with some immune checkpoint genes involved in PD1/CTLA4 pathways. The immune adhesion molecule CD58, also termed lymphocyte function-associated antigen-3 (LFA-3), is a costimulatory receptor extensively expressed on human cells [54, 55]. The interaction between CD58 and its natural ligand (CD2) promotes optimal T/nature killer cell activation and triggers a series of intracellular signaling. Accumulating evidence has demonstrated the central role of CD2-CD58 interaction in modulating antiviral responses, inflammatory responses, and immune evasion of solid cancer cells [55, 56]. In gastric cancer, elevated expression of CD58 was associated with deteriorated tumor cell invasion, reduced survival time, and cancer recurrence [57]. Consistently, our study showed that CD58 might be an unfavorable prognostic gene, a higher expression of which was significantly correlated with shorter disease-free and overall survival. We found that CD58 was highly expressed in tumor tissues with high infiltrating levels of immune cells, suggesting that the activated infiltrating immune cells in tumor tissues might also be suppressed by CD58. Interestingly, CD58 loss was observed to induce immune evasion in multiple melanoma cell and tumor-infiltrating lymphocyte coculture models, and the expression of CD58 was downregulated in melanoma patients receiving immune checkpoint inhibitors [56]. Also, PD-L1 was upregulated in CD58-knockout melanoma cells. Their results suggested that CD58 downregulation (or loss) might contribute to immune evasion via multiple distinct mechanisms (e.g., increased expression of coinhibitory PD-L1) [56]. Additional studies are necessary to investigate the role of CD58 in HCC and CD58-PD-L1 balance.

Furthermore, large amounts of lncRNAs have been demonstrated to be abnormally expressed in HCC, and they play essential roles in cancer development, proliferation, and differentiation [58-60]. LncRNAs could regulate the expression of mRNA in various ways both in cis, and trans [61], and lncRNA can also act as sponges of mRNA or work with protein to regulate the expression of mRNA [62, 63]. Our results shed light on the interaction of mRNA and lncRNA and indicated the important roles of lncRNAs in HCC. Still, the lack of confirmatory experiments is a significant limitation, and more studies are required to validate the results further.

5. Conclusion

Our research identified key functional modules and immunomodulatory regulators for HCC, which might offer novel diagnostic biomarkers and/or therapeutic targets for cancer immunotherapy.