Skip to main content

Screening and analysis of RNAs associated with activated memory CD4 and CD8 T cells in liver cancer

Abstract

Background

Liver cancer is one of the most common malignant tumors in the world. T cell-mediated antitumor immune response is the basis of liver cancer immunotherapy.

Objective

To screen and analyze the RNAs associated with activated memory CD4 T cells and CD8 T cells in liver cancer.

Methods

ESTIMATE was used to calculate the stromal and immune scores of tumor samples, which were downloaded from The Cancer Genome Atlas (TCGA). The differentially expressed genes (DEGs) in high and low stromal and immune scores were screened, followed by functional enrichment of overlapped DEGs. We then conducted a survival analysis to identify immune-related prognostic indicators and constructed protein-protein interaction (PPI) networks and ceRNA networks. Finally, chemical small-moleculetarget interaction pairs associated with liver cancer were screened.

Results

A total of 55,955 stromal-related DEGs and 1811 immune-related DEGs were obtained. The 1238 overlapped DEGs were enriched in 1457 biological process terms and 74 KEGG pathways. In addition, a total of 120 activated memory CD4 T cell-related genes and 309 CD8 T cell-related genes were identified. The survival analysis revealed that upregulated expression of T cell-related genes including EOMES, CST7, and CD5L indicated the favorable prognosis of liver cancer. EOMES was regulated by has-miR-23b-3p and has-miR-23b-3p was regulated by lncRNA AC104820.2 in the ceRNA network of activated memory CD4 T cell-related genes. In addition, EOMES was regulated by has-miR-23a-3p and has-miR-23a-3p was regulated by lncRNA AC000476.1 in the ceRNA network of CD8 T cells.

Conclusion

T cell-related RNAs EOMES, CST7, CD5L, has-miR-23b-3p, and has-miR-23a-3p may be associated with the prognosis of liver cancer. And the molecular characteristics of these T cell-related genes were plotted.

Highlights

A total of 309 CD8 T cell-related genes and 120 activated memory CD4 T cell-related genes were screened in liver cancer tumor samples.

Forty-four chemical small-molecule–target interaction pairs associated with activated memory CD4 T cells and 276 pairs associated with CD8 T cells were screened.

Upregulated expression of T cell-related genes including EOMES, CST7, and CD5L indicated the favorable prognosis of liver cancer.

Introduction

Liver cancer is one of the ordinary malignant tumors, with high morbidity and mortality [1, 2]. Although the incidence of liver cancer has decreased, there are still approximately 840,000 new cases worldwide every year [3]. Surgical resection and liver transplantation are currently recognized as effective methods for the treatment of liver cancer; however, the 5-year survival rate is only approximately 50%, and the 5-year recurrence rate can reach 60~70% [4]. Epigenetic-related molecules (such as SNRPC) [5], non-coding RNA (such as hsa-mir-221) [6], and immune-related molecules (such as DCK) [7] could be used as potential biomarkers for diagnosis, treatment, and prognosis monitoring in liver cancer. Therefore, it is critically necessary to further explore the molecular pathogenesis of this disease by screening for new molecular biomarkers that could enhance early diagnosis and improve current therapies for liver cancer.

According to the competing endogenous RNA (ceRNA) hypothesis, different RNA transcripts participate in pathological processes, primarily by competing for the binding sites of shared miRNAs. Li et al. have suggested that long noncoding RNA LINC01139 promotes the progression of liver cancer by upregulating MYBL2 via competitively binding to members of the miR-30 family [8]. In addition, tumor infiltration of immune cells is closely related to the prognosis of tumors and the determination of immunotherapeutic targets [9]. Tumor immune escape is an important characteristic of tumor formation and is related to the decline of the T cell response-ability [10]. T cells are key mediators of tumor destruction and are critical to the specificity of tumor antigen expression [11]. Moreover, numerous studies have shown that the T cell-mediated antitumor immune response is the basis of tumor immunotherapy, which is associated with a favorable prognosis [12]. Parriott et al. found that T cells expressing the chimeric-PD1-Dap10-CD3zeta receptor can reduce the tumor burden in a variety of mouse solid cancer syngeneic models [13]. CD8+ T lymphocytes are essential factors affecting the efficacy of immunotherapy. Donghua et al. identified CD8+ T cell co-expression genes (C1QC, CD3D, GZMA, and PSMB9 ) that promoted infiltration of CD8+ T cells in liver cancer [14]. Immunotherapy is widely used to treat advanced cancer. However, immunotherapy has many limitations, such as low success rate, many complications, and rapid progression. Studies have shown that the low success rate of immunotherapy is related to less T lymphocyte infiltration. Analysis of infiltrating T cell related molecules may be helpful for immunotherapy.

However, there are few reports on the immune-related prognostic indicators of liver cancer. In this study, we aimed to analyze T cells related to tumor immunity and investigate RNAs associated with T cells in liver cancer to identify useful molecular markers related to the prognosis of liver cancer. In the tumor microenvironment, immune cells and stromal cells are two main types of non-tumor components. They are of great significance for the diagnosis and prognosis evaluation. The immune and stroma scores which was based on the ESTIMATE algorithm are the quantification of the immune and stroma components in the tumor microenvironment. The differential genes were divided into upregulated group and downregulated group according to the intermediate value of quantitative scores. In order to find more representative genes, common stroma and immune genes were screened in the upregulated group and downregulated group, respectively. Then, we analyzed the molecular mechanisms associated with T cells in liver cancer based on the representative genes.

Materials and methods

Data sources and data preprocessing

The RNA expression RNAseq sequencing data (Counts) and clinical data (phenotype) of LIHC (Version: 07-19-2019) were downloaded from the TCGA database [15] (https://xenabrowser.net/). Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) was used to standardize the original data. The FPKM of a gene in sample was equal to the ratio of the total exon fragments to the reading falling in the genome (map reading (Millions)) and the length of the gene (exon length (KB)). FPKM = Total exon fragments/(Mapped reads×Exon length). The RNA expression RNAseq sequencing data (FPKM) was also downloaded from the TCGA database. A total of 58,387*369 Counts and FPKM representation matrices were obtained. The detection platform of sequencing data was the Illumina HiSeq 2 000 RNA Sequencing platform.

The RNA-Seq was annotated based on the annotation file of the Gencode database [16] (https://www.gencodegenes.org/). Then, genes with “protein_coding” annotations were extracted as mRNA. Meanwhile, RNAs with “TEC,” “known_ncrna,” “macro_lncRNA,” “bidirectional_promoter_lncrna,” and “lncRNA” annotations were extracted as lncRNA, and other genes were defined as “undefined.” Then, the sequencing data log2 (count+1) values of Counts were restored to the count value. Reads and the sequencing data log2 (count+1) values of FPKM were restored to the fpkm value. In the expression profile analysis of the Ensembl_ID, the mapping probe was used to calculate the gene expression value (obtained from the annotation files of the chip platform and microarray dataset) to Symbol_ID. The average value was taken as the level of Ensembl_ID expression when multiple probes matched one Symbol_ID. Lastly, 369 tumor samples, 60,483 ENSEMBLE IDs, and 58,387 RNA symbols were obtained in Table 1. The clinical information related to prognosis in TCGA was also collected in Table 1, containing overall survival (OS) and OS status, and the lifetime was converted from days to months (days/30).

Table 1 The basic statistical information of tumor samples in the TCGA-LIHC dataset

Immune score calculation

The stromal and immune scores of all tumor samples were calculated with the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE algorithm in R [17].

Differential expression analysis

The tumor samples were split into high- and low-score groups based on the median value of stromal and immune scores of all tumor samples. The typical Bayesian modified t test in limma package [18] (Version 3.40.6) was applied to analyze differentially expressed RNA (DERNA) between high and low immune-score groups with a cutoff of |logFC| > 1 and P < 0.05. Lastly, the ggscatter function of the ggpubr package [19] in the R language (Version: 0.2.2) was used to draw a volcano plot of the DERNAs. The DEGs of stromal and immune-score groups were selected as the DEGs related to the immune microenvironment.

Enrichment analysis

Based on up- and downregulated DEGs, the clusterProfiler package [20] in R was used to perform Gene Ontology (GO) [21] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [22] pathway enrichment analysis with a cutoff of P < 0.05 and count ≥ 2.

Screening of T cell-related genes

RNA-Seq expression profile data were used to target DEGs, and the abundance matrix of the immune cells of the samples was evaluated using the CIBERSORT deconvolution algorithm [23] to analyze the abundance of infiltration of 22 immune cells in tumor samples. The LM22 dataset provided by the Cibersort website was used as the RNA expression signature template following the analysis of the abundance of infiltration of immune cells, and parameters were set as perm = 50 and QN = N. The landscape of the immune cells was then plotted using the ggplot2 package in R. The Pearson correlation coefficient between the expression value of the DEGs and the abundance of immune infiltration in activated memory CD4 T cells and CD8 T cells was calculated, and the activated memory CD4 T and CD8 T cell-related genes were corrected for BH with a threshold value of |r| > 0.3 and P < 0.05. Moreover, KEGG pathway analysis on T cell-related genes was carried out based on the carcinoma- and hepatocellular-related pathways in the Comparative Toxicogenomics Database (CTD) [24] (http://ctdbase.org/).

Survival analysis

Based on the median expression in T cell-related genes, samples were divided into high- and low-expression groups. The survival analysis was then carried out using a Kaplan-Meier plot. Finally, a log-rank statistical test was conducted with a significance threshold of P < 0.05.

PPI network analysis

The STRING database [25] was used to analyze the protein-protein interactions encoded by T cells. The PPI score was set as 0.7 (high-confidence value). Afterward, the PPI network of T cell-related genes was constructed using Cytoscape software [26].

Prediction of miRNA-target interaction pairs

The miRNAs of the T cell-related genes were predicted using miRWalk 3.0 [27]. Then, the miRNA-target interaction pairs in the TargetScan [28], MiRDB [29], and MirTarBase [30] databases were also obtained using a threshold of score > 0.95. In addition, the HMDD V3.2 database [31] was used to further validate and screen the predicted miRNA with the keywords “Carcinoma, Hepatocellular”.

Construction of the ceRNA network

The lncRNAs-miRNAs relationship between T cell-related genes was predicted using the DIANA-LncBase database [32]. Finally, the T cell-related gene lncRNA-miRNA-mRNA network was constructed using Cytoscape software.

Chemical small-molecule–target network analysis

To search for liver cancer-related genes and chemicals, the Comparative Toxicogenomics Database (http://ctd.mdibl.org/) was searched using “Carcinoma, Hepatocellular” as keywords. Genes that were both associated with liver cancer and belonged to the T cell-related genes ceRNA network were used to screen chemical-target pairs. The T cell-related genes chemical small-moleculetarget network was obtained utilizing the Cytoscape software.

Results

The DEGs between the high and low stromal and immune score groups

The aim was to screen immune-related differential genes in liver cancer. In total, 10,221 genes from the TCGA-LIHC dataset were matched by the ESTIMATE algorithm, and the number of unmatched genes was 191. There were 1811 DEGs (of which 1744 and 67 were up- and downregulated, respectively) were selected according to the high and low immune score groups (Fig. 1A and Supplementary file 1). Meanwhile, a total of 55,955 DEGs (of which 2279 and 153 were up- and downregulated, respectively) were selected between high and low stromal score groups (Fig. 1B and Supplementary file 2). In addition, 1211 overlapped upregulated DEGs and 27 overlapped downregulated DEGs were selected in the stromal score and immune score groups (Fig. 1C, D). The names of DEGs are shown in Supplementary file 3. In the enrichment analysis of 1238 overlapped DEGs, the GO analysis showed that the 1238 DEGs were mainly enriched in 1457 GO-biological processes (BPs; e.g., regulation of lymphocyte activation, leukocyte migration, and lymphocyte differentiation) and 74 KEGG pathways [e.g., cell-adhesion molecules (CAMs)] (Fig. 1E, F).

Fig. 1
figure 1

Differentially expressed genes (DEGs) between the high and low stromal and immune score groups. Based on the median value of stromal and immune scores of all tumor samples, the tumor samples were divided into high- and low-score groups. The standard Bayesian modified t test in the limma package was applied to analyze differentially expressed RNAs (DERNA) between high and low immune-score groups with a cutoff of |logFC| > 1 and P < 0.05. According to the stromal score, there were 55,955 DEGs (2279 upregulated and 153 downregulated) were selected between high and low stromal-score groups. Meanwhile, 1811 DEGs (including 1744 upregulated and 67 downregulated) were selected according to the high and low immune-score groups. In addition, there were 1211 upregulated overlapped DEGs and 27 downregulated overlapped DEGs were selected in the stromal score and immune score groups. In the enrichment analysis of 1238 overlapped DEGs, the GO analysis showed that the 1238 DEGs were mainly enriched in 1457 GO-BPs and 74 KEGG pathways. A The DEGs between high and low stromal-score groups. B The DEGs between high and low immune-score groups. C The upregulated DEGs in the stromal score and immune score groups. D: The downregulated DEGs in the stromal score and immune score groups. E The GO analysis of all overlapped DEGs. F The KEGG pathways analysis of all overlapped DEGs. Red nodes represent upregulated DEGs and blue nodes represent downregulated DEGs. The size of the ball represents the number of genes enriched in each term. The color of the ball represents the value of the P value. DEGs: differentially expressed genes; GO, Gene Ontology; BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes

Identification of T cell-related genes

T cells played an important role in liver cancer occurrence and development. The aim was to further analyze T cell-related genes in liver cancer. The abundance of immune-cell infiltration in tumor samples was estimated and the immune-cell landscape was shown in Fig. 2. Meanwhile, according to the Pearson correlation coefficient, a total of 120 activated memory CD4 T cell-related genes (Supplementary file 4) and 309 CD8 T cell-related genes were identified (Supplementary file 5). In the enrichment analysis of activated memory CD4 T cells and CD8 T cell-related genes, activated memory CD4 T cell-related genes were enriched in 406 GO-BPs and 35 KEGG pathways [e.g., CAMs, Th1 and Th2 cell differentiation, and Hematopoietic cell lineage], the Top 8 GO terms and Top 10 KEGG pathways were shown in Fig. 3A and B, respectively. In addition, CD8 T cell-related genes were involved in 596 GO-BPs and 42 KEGG pathways [e.g., Th17 cell differentiation, Th1 and Th2 cell differentiation, and CAMs], the Top 8 GO terms and Top 10 KEGG pathways were shown in Fig. 3C and D, respectively. In addition, based on the CTD database, activated memory CD4 T cell-related genes were involved in 31 liver cancer-related pathways [e.g., CAMs, Th1 and Th2 cell differentiation, hematopoietic cell lineage] (Table 2), and CD8 T cell-related genes were involved in 36 liver cancer-related pathways [e.g., Th17 cell differentiation, Th1 and Th2 cell differentiation, and CAMs] (Table 3).

Fig. 2
figure 2

The landscape of the immune cells. To analyze the abundance of infiltration of the immune cells in the samples, RNA-Seq expression profile data were used to target DEGs, and the abundance matrix of the immune cells was evaluated through the CIBERSORT deconvolution algorithm. Finally, the abundance of infiltration of immune cells (including naïve B cells, memory B cells, CD8 T cells, activated memory CD4T cells, M0 macrophages, M1 macrophages, M2 macrophages, activated dendritic cells, and neutrophils) was determined

Fig. 3
figure 3

Enrichment analysis of activated memory CD4 T cells and CD8 T cell-related genes. The clusterProfiler package in R was used to perform GO and KEGG pathway enrichment analysis with a cutoff of P < 0.05 and count ≥ 2. In the enrichment analysis of activated memory CD4 T cells and CD8 T cell-related genes, activated memory CD4 T cell-related genes were enriched in 406 GO-BPs (e.g., regulation of lymphocyte activation, regulation of T cell activation, and regulation of cell-cell adhesion) and 35 KEGG pathways [e.g., CAMs, Th1 and Th2 cell differentiation and Hematopoietic cell lineage]. In addition, CD8 T cell-related genes were involved in 596 GO-BPs (e.g., regulation of lymphocyte activation, regulation of T cell activation, and regulation of cell-cell adhesion) and 42 KEGG pathways [e.g., Th17 cell differentiation, Th1 and Th2 cell differentiation, and CAMs]. A The GO analysis of activated memory CD4 T cell-related genes. B The GO analysis of CD8 T cell-related genes. C The KEGG pathways analysis activated memory CD4 T cell-related genes. D The KEGG pathways analysis CD8 T cell-related genes. The size of the ball represents the number of genes enriched in each term. The color of the ball represents the value of the P value. GO, Gene Ontology; BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes

Table 2 The KEGG pathway analysis of activated memory CD4 T cells and CD8 T cell-related genes in the Comparative Toxicogenomics Database. The KEGG pathway analysis of activated memory CD4 T cell-related genes
Table 3 The KEGG pathway analysis of activated memory CD4 T cells and CD8 T cell-related genes in the Comparative Toxicogenomics Database. The KEGG pathway analysis of CD8 T cell-related genes

T cells associated with prognosis of liver cancer

Liver cancer was an often fatal malignant with poor prognosis. The purpose was to further analyze the role of T cell-related genes in the prognosis of hepatocellular carcinoma. A total of 363 sets of survival information from patients with liver cancer was obtained. Moreover, nine activated memory CD4 T cell-related genes were associated with liver cancer prognosis [e.g., eomesodermin (EOMES), glutathione S-transferase (CST7), and adhesion G protein-coupled receptor E2 (EMR2)] (Table 4). In addition, there were 30 CD8 T cell-related genes associated with liver cancer prognosis [e.g., CD5 molecules like CD5L, EOMES, and CST7] (Table 5). Upregulated expression of T cell-related genes including EOMES, CST7, and CD5L indicated the favorable prognosis of liver cancer. Downregulated expression of T cell-related genes including NCF2 and HTAR3 indicated the poor prognosis of liver cancer.

Table 4 Activated memory CD4 T cells and CD8 T cells associated with prognosis of liver cancer. Activated memory CD4 T cell-related genes associated with liver cancer prognosis
Table 5 Activated memory CD4 T cells and CD8 T cells associated with prognosis of liver cancer. CD8 T cell-related genes associated with liver cancer prognosis

PPI network of the T cell-related genes

The PPI network was constructed to plot the characteristics of these molecules. PPI network analysis of activated memory CD4 T cell-related genes revealed 53 nodes, including one survival-related gene (EOMES), and 162 interaction pairs (Fig. 4A). The CD8 T cell-related gene PPI network contained 127 nodes, including 11 survival-related genes [e.g., EOMES, CD69 molecule (CD69), and zeta chain of T cell receptor-associated protein kinase 70 (ZAP70)], and 613 interaction pairs (Fig. 4B).

Fig. 4
figure 4

The PPI network of activated memory CD4 and CD8 T cell-related genes. The STRING database was used to analyze protein-protein interactions encoded by activated memory CD4 T cells and CD8 T cells. The PPI score was set as 0.7 (high-confidence value). Afterward, the PPI networks of activated memory CD4 T cells and CD8 T cell-related genes were constructed using Cytoscape software. PPI network analysis of activated memory CD4 T cell-related genes revealed 53 nodes and 162 interaction pairs. The CD8 T cell-related gene PPI network contained 127 nodes and 613 interaction pairs. A The PPI network of activated memory CD4 T cell-related genes. B The PPI network of CD8 T cell-related genes. Red nodes represent survival-related DEGs, triangles represent upregulated DEGs, and blue nodes represent other DEGs. The size of nodes represents the value. Larger nodes indicate a larger value. PPI, protein-protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes

CeRNA network of the T cell-related genes

The ceRNA network was constructed to describe the regulatory effect of non-coding RNA on T cell-related differential molecules. Based on the HMDD database, ten miRNA-mRNA relationships of activated memory CD4 T cell-related genes were obtained (10 miRNAs and three target genes). Also, 26 miRNA-mRNA relationships of CD8 T cell-related genes were obtained (22 miRNAs and 10 target genes). Four miRNA-lncRNA relationships of activated memory CD4 T cell-related genes were obtained (four miRNAs and one lncRNA). Furthermore, 21 miRNA-lncRNA relationships of CD8 T cell-related genes were obtained (21 miRNAs and 13 lncRNAs). Based on the ten miRNA-mRNA relationships and four miRNA-lncRNA relationships of activated memory CD4 T cell-related genes, a ceRNA network of activated memory CD4 T cell-related genes was constructed (Fig. 5A). Here, EOMES was regulated by has-miR-23b-3p and has-miR-23b-3p was regulated by lncRNA AC104820.2. In addition, the ceRNA network of CD8 T cell-related genes was constructed among 26 miRNA-mRNA relationships and 21 miRNA-lncRNA relationships (Fig. 5B). Here, EOMES was regulated by has-miR-23a-3p and has-miR-23a-3p was regulated by lncRNA AC000476.1.

Fig. 5
figure 5

CeRNA network of activated memory CD4 T cells and CD8 T cell-related genes. The miRNAs of activated memory CD4 T cells and CD8 T cell-related genes were predicted using miRWalk 3.0, and miRNA-target interaction pairs in the TargetScan, MiRDB, and MirTarBase databases were obtained using a threshold of score > 0.95. In addition, the HMDD V3.2 database was used to further validate and screen the predicted miRNA using the keywords “Carcinoma, Hepatocellular”. The lncRNAs-miRNAs relationship between activated memory CD4 T cells and CD8 T cell-related genes was predicted using the DIANA-LncBase database. Finally, activated memory CD4 T cell- and CD8 T cell-related gene lncRNA-miRNA-mRNA network was constructed utilizing Cytoscape software. EOMES was regulated by has-miR-23b-3p and has-miR-23b-3p were regulated by lncRNA AC104820.2. EOMES was regulated by has-miR-23a-3p and has-miR-23a-3p was regulated by lncRNA AC000476.1. A ceRNA network of activated memory CD4 T cell-related genes. B ceRNA network of CD8 T cell-related genes. Red nodes represent upregulated DEGs, green triangles represent miRNA, and the red rhombi represent upregulated lncRNAs

Chemical small-molecule–target network analysis of T cell-related genes

The chemical small-molecule–target network was constructed to find agents that regulate these differentially expressed genes. There were 44 chemical small-moleculetarget interaction pairs associated with activated memory CD4 T cells (Fig. 6A), including five mRNAs and 26 chemical small molecules. In addition, there were 276 CD8 T cell-associated chemical small-moleculetarget interaction pairs, containing 19 mRNAs and 110 chemical small molecules (Fig. 6B).

Fig. 6
figure 6

Chemical small-molecule–target network analysis of activated memory CD4 T cells and CD8 T cell-related genes. To search for liver cancer-related genes and chemicals, the Comparative Toxicogenomics Database was searched using “Carcinoma, Hepatocellular” as keywords. Genes that were both associated with liver cancer, and belonged to the T cell-related genes ceRNA network, were used to screen chemical-target pairs. The T cell-related genes chemical small-moleculetarget network was obtained utilizing the Cytoscape software. There were 44 chemical small-moleculetarget interaction pairs associated with activated memory CD4 T cells, including five mRNAs and 26 chemical small molecules. In addition, there were 276 CD8 T cell-associated chemical small-moleculetarget interaction pairs, containing 19 mRNAs and 110 chemical small molecules. A The chemical small-moleculetarget network of genes in activated memory CD4 T cells. B The chemical small-moleculetarget network of genes in CD8 T cells. Red nodes represent survival-related upregulated DEGs, and green nodes represent chemical small molecules. The size of nodes represents the value, such that larger nodes indicate a larger value

Discussion

In this study, a total of 55,955 stromal-related DEGs and 1811 immune-related DEGs were obtained. Then, the 1238 overlapped DEGs were enriched in 1457 BPs and 74 KEGG pathways. In addition, a total of 120 activated memory CD4 T cell-related genes and 309 CD8 T cell-related genes were identified. The survival analysis showed that the T cell-related genes EOMES, CST7, CD5L, and EMR2 were associated with the prognosis of liver cancer. Our analysis recommended the AC104820.2-has-miR-23b-3p-EOMES axis and the AC000476.1-has-miR-23a-3p-EOMES axis from activated memory CD4 T cell- and CD8 T cell-related ceRNAs, respectively.

A total of 120 activated memory CD4 T cell-related genes and 309 CD8 T cell-related genes were identified in this study. Moreover, T cell-related genes were involved in Th1 and Th2 cell differentiation and CAM signaling pathways. Wang et al. found that abnormal expression of CD4 T cell subsets in peripheral blood of patients with ovarian cancer and an imbalance of Th1/Th2 and Treg/Th17 provide a basis for clinical immunotherapy of ovarian cancer [33]. In addition, Su et al. identified lncRNAs and genes related to the prognosis of gastric cancer and showed that prognostic genes positively related to mortality were enriched in the CAM signaling pathway by performing a functional enrichment analysis [34]. Furthermore, Talia et al. showed that GPER levels are related to metastasis and premigration genes belonging to the CAM signaling pathway in breast cancer [35], which further suggested that T cell-related genes may be related to the progression of liver cancer via Th1- and Th2-cell differentiation and CAM signaling pathways.

EOMES is an important transcription factor that regulates the differentiation and function of effector T cells, and previous studies have shown that the survival of tumor patients is closely related to EOMES expression. For instance, Gao et al. suggested that the potential anticancer functions of EOMES, ATF5, and ECM1 were confirmed by siRNA experiments [36]. The CD8 T cells were divided into three distinct subpopulations: PD1Hi, PD1Int, and PD1-. Ma et al. showed that compared with adjacent non-tumor liver tissues, the PD1Hi CD8 T cells were significantly enriched in tumors, and the PD1Hi CD8 T cells in liver cancer were highly expressed with depletion-related inhibitory receptors (TIM3, ctla-4, and so on) and transcription factors (EOMES, BATF, and so on) [37]. In the present study, EOMES was a T cell-related gene associated with liver cancer prognosis. Moreover, EOMES was regulated by has-miR-23b-3p and has-miR-23b-3p was regulated by lncRNA AC104820.2 in activated memory CD4 T cell-related genes ceRNA network, and EOMES was regulated by has-miR-23a-3p and has-miR-23a-3p was regulated by lncRNA AC000476.1 in CD8 T cell-related genes ceRNA network. Zaman et al. found an inhibitory effect of miR-23b-3p on the expression of the PTEN gene in renal carcinoma [38]. Meanwhile, Wen et al. showed that osthole inhibited EMT-mediated metastasis of prostate cancer by inhibiting snail signaling and miR-23a-3p [39]. Considering these results together, we speculate that EOMES may be the potential target of has-miR-23b-3p and has-miR-23a-3p in liver cancer and that they play important roles in the progression of the disease.

Moreover, in the present study, CST7 was a T cell-related gene associated with liver cancer prognosis. CST7, also known as CMAP, has been reported to have close connections with liver cancer. For instance, Zhou et al. suggested that CST7 and CSTB genes may serve as potential prognostic and diagnostic biomarkers for liver cancer [40]. Interestingly, CMAP is a novel cystatin-like gene involved in liver metastasis [41]. CD5L is a soluble scavenger cysteine-rich protein that regulates the inflammatory response. Aran et al. suggest that CD5L is upregulated in hepatocellular carcinoma and promotes the proliferation and anti-apoptotic response of liver cancer cells by binding to HSPA5 (GRP78) [42]. Thus, we speculate that CST7 and CD5L contribute to liver cancer progression.

In the present study, immune-related prognostic indicators and chemical small-moleculetarget interaction pairs associated with liver cancer were screened, which could guide liver cancer clinical decision-making. Moreover, we constructed a ceRNA network of immune-related genes associated with liver cancer prognosis and predicted miRNA targets, which provided a basis for the study of tumor-induced T cell-mediated immune escape. However, we did not use molecular biological experiments to verify our results. The molecular mechanisms that T cell-related genes influence the prognosis of liver cancer have not been investigated, and further molecular biological experiments are needed to explore the specific roles of these genes.

Conclusion

T cell-related RNAs EOMES, CST7, CD5L, has-miR-23b-3p, and has-miR-23a-3p may be associated with the prognosis of liver cancer. And the molecular characteristics of these T cell-related genes were plotted by PPI network, ceRNA network, chemical small-molecule–target network.

Availability of data and materials

Not applicable.

Abbreviations

TCGA:

The Cancer Genome Atlas

KEGG:

Kyoto Encyclopedia of Genes and Genomes

GO:

Gene ontology

PPI:

Protein-protein interaction

BP:

Biological process

ESTIMATE:

Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data

CD5L :

CD5 molecule like

EOMES :

Eomesodermin

CST7 :

Glutathione S-transferase

CD69 :

CD69 molecule

ZAP70 :

Zeta chain of T cell receptor-associated protein kinase 70

References

  1. Lee S, et al. Diagnostic performance of CT/MRI liver imaging reporting and data system v2017 for hepatocellular carcinoma: a systematic review and meta-analysis. Liver Int. 2020;40(6):1488–97.

  2. Zhiyong D. Cinobufacini injection for moderate and advanced primary liver cancer: a systematic review and meta-analysis. J Chin Pharm Sci. 2019;28(4):264–75.

    Article  Google Scholar 

  3. Valery PC, et al. Projections of primary liver cancer to 2030 in 30 countries worldwide. Hepatology. 2018;67(2):600–11.

    Article  Google Scholar 

  4. Margonis GA, et al. Association of BRAF mutations with survival and recurrence in surgically treated patients with metastatic colorectal liver cancer. JAMA Surg. 2018;153(7):e180996.

    Article  Google Scholar 

  5. Jianxin X, et al. N6-methyladenosine (m6A) RNA methylation regulator SNRPC is a prognostic biomarker and is correlated with immunotherapy in hepatocellular carcinoma. World J Surg Oncol. 2021;19(1):241.

    Article  Google Scholar 

  6. Zhongjun W, et al. Identification and interaction analysis of key genes and microRNAs in hepatocellular carcinoma by bioinformatics analysis. World J Surg Oncol. 2017;15(1):63.

    Article  Google Scholar 

  7. Wang X, et al. DCK is a promising prognostic biomarker and correlated with immune infiltrates in hepatocellular carcinoma. World J Surg Oncol. 2020;18(1):176.

    Article  Google Scholar 

  8. Li ZB, et al. Long noncoding RNA LINC01139 promotes the progression of hepatocellular carcinoma by upregulating MYBL2 via competitively binding to miR-30 family. Biochem Biophys Res Commun. 2020;525(3):581–8.

  9. Rohr-Udilova N, et al. Deviations of the immune cell landscape between healthy liver and hepatocellular carcinoma. Sci Rep. 2018;8(1):6220.

    Article  CAS  Google Scholar 

  10. Dmitry G. Mechanisms and functional significance of tumour-induced dendritic-cell defects. Nat Rev Immunol. 2004;4(12):941–52.

    Article  CAS  Google Scholar 

  11. Kishton RJ, Sukumar M, Restifo NP. Metabolic regulation of T cell longevity and function in tumor immunotherapy. Cell Metab. 2017;26(1):94–109.

    Article  CAS  Google Scholar 

  12. Ribas A, Wolchok JD. Cancer immunotherapy using checkpoint blockade. Science. 2018;359(6382):1350–5.

  13. Parriott G, et al. T cells expressing a chimeric-PD1-Dap10-CD3zeta receptor reduce tumor burden in multiple murine syngeneic models of solid cancer. Immunology. 2020;160(3):280–94.

  14. Qi Pan Y, et al. Identification of CD8+ T cell-related genes: correlations with immune phenotypes and outcomes of liver cancer. J Immunol Res. 2021;2021:9960905.

    PubMed  PubMed Central  Google Scholar 

  15. Tomczak K, Czerwińska P, Wiznerowicz M. The cancer genome atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol. 2015;19(1A):A68.

    Google Scholar 

  16. Harrow J, et al. GENCODE: the reference human genome annotation for the ENCODE project. Genome Res. 2012;22(9):1760–74.

    Article  CAS  Google Scholar 

  17. Yoshihara K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4(1):1–11.

    Article  CAS  Google Scholar 

  18. Smyth GK, et al. LIMMA: linear models for microarray data. In: Bioinformatics and computational biology solutions using R and bioconductor. Statistics for biology and health; 2005.

    Google Scholar 

  19. Kassambara A. Ggpubr:“ggplot2” based publication ready plots. R package version 0.1; 2017. p. 6.

    Google Scholar 

  20. Yu G, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.

    Article  CAS  Google Scholar 

  21. Ashburner M, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9.

    Article  CAS  Google Scholar 

  22. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  Google Scholar 

  23. Newman AM, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.

    Article  CAS  Google Scholar 

  24. Davis AP, et al. The comparative toxicogenomics database: update 2019. Nucleic Acids Res. 2019;47(D1):D948–54.

    Article  CAS  Google Scholar 

  25. Szklarczyk D, et al. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res. 2016;45(D1):D362–8.

  26. Shannon P, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  Google Scholar 

  27. Dweep H, et al. miRWalk–database: prediction of possible miRNA binding sites by “walking” the genes of three genomes. J Biomed Inform. 2011;44(5):839–47.

    Article  CAS  Google Scholar 

  28. Agarwal V, et al. Predicting effective microRNA target sites in mammalian mRNAs. elife. 2015;4:e05005.

    Article  Google Scholar 

  29. Wong N, Wang X. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res. 2015;43(D1):D146–52.

    Article  CAS  Google Scholar 

  30. Chou C-H, et al. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018;46(D1):D296–302.

    Article  CAS  Google Scholar 

  31. Huang Z, et al. HMDD v3. 0: a database for experimentally supported human microRNA–disease associations. Nucleic Acids Res. 2019;47(D1):D1013–7.

    Article  CAS  Google Scholar 

  32. Paraskevopoulou MD, et al. DIANA-LncBase v2: indexing microRNA targets on non-coding transcripts. Nucleic Acids Res. 2016;44(D1):D231–8.

    Article  CAS  Google Scholar 

  33. Wang LH, et al. Th1/Th2 and Treg/Th17 cell balance in peripheral blood of patients with ovarian cancer. Nan Fang Yi Ke Da Xue Xue Bao. 2017;37(8):1066–70.

    PubMed  CAS  Google Scholar 

  34. Su X, et al. Identification of the prognosis-related lncRNAs and genes in gastric cancer. Front Genet. 2020;11:27.

    Article  CAS  Google Scholar 

  35. Talia M, et al. The G protein-coupled estrogen receptor (GPER) expression correlates with pro-metastatic pathways in ER-negative breast cancer: a bioinformatics analysis. Cells. 2020;9(3):622.

  36. Gao F, et al. Integrated analyses of DNA methylation and hydroxymethylation reveal tumor suppressive roles of ECM1, ATF5, and EOMES in human hepatocellular carcinoma. Genome Biol. 2014;15(12):533.

    Article  CAS  Google Scholar 

  37. Ma J, et al. PD1(hi) CD8(+) T cells correlate with exhausted signature and poor clinical outcome in hepatocellular carcinoma. J Immunother Cancer. 2019;7(1):331.

    Article  Google Scholar 

  38. Zaman MS, et al. Inhibition of PTEN gene expression by oncogenic miR-23b-3p in renal cancer. PLoS One. 2012;7(11):e50203.

    Article  CAS  Google Scholar 

  39. Wen YC, et al. By inhibiting snail signaling and miR-23a-3p, osthole suppresses the EMT-mediated metastatic ability in prostate cancer. Oncotarget. 2015;6(25):21120–36.

    Article  Google Scholar 

  40. Zhou X, et al. Investigation of the clinical significance and prospective molecular mechanisms of cystatin genes in patients with hepatitis B virusrelated hepatocellular carcinoma. Oncol Rep. 2019;42(1):189–201.

    PubMed  PubMed Central  CAS  Google Scholar 

  41. Morita M, et al. CMAP: a novel cystatin-like gene involved in liver metastasis. Cancer Res. 1999;59(1):151–8.

    PubMed  CAS  Google Scholar 

  42. Aran G, et al. CD5L is upregulated in hepatocellular carcinoma and promotes liver cancer cell proliferation and antiapoptotic responses by binding to HSPA5 (GRP78). FASEB J. 2018;32(7):3878–91.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge the multiple databases, which made the data available.

Funding

This work was supported by the Welfare Technology Application Research Program of Huzhou (No. 2019GY01).

Author information

Authors and Affiliations

Authors

Contributions

All authors participated in the conception and design of the study; conceived and drafted the manuscript: Zhang Yan, Yin Lijuan, and Wu Yinhang; collated and proofread the literature: Zhang Yan, Pan Yuefen, Wu Wei, and Xu Jiamin; wrote the paper: Jin Yin and Han Shuwen; all authors read and approved the paper.

Corresponding authors

Correspondence to Pan Yuefen or Han Shuwen.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Yan, Z., Lijuan, Y., Yinhang, W. et al. Screening and analysis of RNAs associated with activated memory CD4 and CD8 T cells in liver cancer. World J Surg Onc 20, 2 (2022). https://doi.org/10.1186/s12957-021-02461-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12957-021-02461-6

Keywords

  • Liver cancer
  • Differentially expressed analysis
  • Immune-related genes
  • T cells
  • Survival analysis
  • CeRNA network