Bioinformatic-based mechanism identification of E2F1-related ceRNA and E2F1 immunoassays in hepatocellular carcinoma.

Background: E2F1 is an important transcription factor. Previous studies have shown that the overexpression of E2F1 is closely related to the occurrence and development of hepatocellular carcinoma (HCC). However, the current research on the regulatory mechanism of E2F1 is still insufficient. This study sought to identify valuable therapeutic E2F1-related targets for HCC.
Methods: HCC-related transcriptome data and patient clinical information downloaded from The Cancer Genome Atlas (TCGA) database. The expression of the E2F1 gene in pan-cancer was analyzed using the Tumor IMmune Estimation Resource (TIMER) 2.0 database, and the expression level of E2F1 in HCC was verified using the Gene Expression Profiling Interactive Analysis database. The overall survival (OS) and progression-free survival (PFS) in HCC patients were also analyzed. Subsequently, based on the Encyclopedia of RNA Interactomes (ENCORI) database, we adopted E2F1 as the research objective and identified the target long non-coding RNAs (lncRNAs) and microRNAs that suggested the competing endogenous RNA (ceRNA) mechanisms related to E2F1. We also performed a correlation analysis of E2F1 using the R language package that contained immune cell and immune checkpoint information. Finally, the drug sensitivity of E2F1 was detected using the R language package, "pRRophetic."
Results: Ultimately, the following 6 potential ceRNA-based pathways targeting E2F1 were identified-lncRNA: LINC01224, PCBP1-AS1, and ITGA9-AS1-miR-29b-3p-E2F1; lncRNA: SNHG7 and THUMPD3-AS1, and LINC02323-miR-29c-3p-E2F1. Cluster of differentiation (CD)4 memory activated T cells, memory B cells, eosinophils, and T follicular helper cells were positively correlated with E2F1 (P<0.05), and monocytes, naïve B cells, and CD4 memory resting T cells were negatively correlated with E2F1 (P<0.05). The immune checkpoint analysis showed that E2F1 was positively correlated with PDCD1, CTLA4, and LAG3 (P>0.2). According to the drug sensitivity analysis, E2F1 may be sensitive to 39 drugs (P<0.05). Conclusions: This study provides a valuable direction for researching transcription factor E2F1, which may be conducive in identifying research targets for HCC-related molecular biological therapy and immunotherapy in future. 2022 Journal of Gastrointestinal Oncology. All rights reserved.

Introduction

Hepatocellular carcinoma (HCC) accounts for about 80% of primary liver cancer cases (1). The major pathogenic causes of HCC include hepatitis [e.g., hepatitis B virus (HBV) and hepatitis C virus (HCV)] alcoholism, smoking, obesity, and congenital inheritance. In the United States, the latest cancer statistics showed that the number of deaths from liver cancer reached 30,230 in 2021 (2). HCV infection is the leading cause of liver cancer in Western countries and causes approximately 1/4 of all the HCC cases. In developing countries (e.g., China), HBV infection is the predominant cause of liver cancer (3). Compared to other cancers, the prognosis of liver cancer is still relatively poor, and it has a 5-year survival rate of only 20% (4). Thus, research urgently needs to be conducted to identify effective biological targets related to liver cancer, especially HCC.

Most long non-coding RNAs (lncRNAs) do not encode proteins (5), and have even been considered junk DNA; however, in-depth research on non-coding RNAs have revealed that many lncRNAs regulate gene expression during or after transcriptional processes. LncRNAs affect a series of pathological and physiological processes by participating in the biological regulation, such as chromosome imprinting, epigenetic regulation, cell proliferation and cell cycle (6,7). Under the recently proposed potential competing endogenous RNA (ceRNA) theory, lncRNA competes to occupy a large number of micro RNAs (miRNAs) in the cell and acts like a sponge to buffer and interfere with the protein encoded by the target gene messenger RNA (mRNA) (8). This kind of mechanism also provides a good entry point for researchers to explore the mechanism of tumorigenesis and development and find effective tumor therapy targets.

E2F1 is an important transcription factor involved in multiple steps, including DNA damage response and cell-cycle regulation (9). Previous studies have shown that the overexpression of E2F1 is closely related to the occurrence and development of various malignant tumors, including HCC (10-12). The abnormal activation of E2F1 affects its downstream transcriptional targets, resulting in DNA replication stress (13). The above mechanisms play an important role in the occurrence and development of liver cancer. At present, there is still a lot of room for exploration on the regulatory mechanism upstream of E2F1. By finding out the effective regulatory mechanism related to E2F1, and then inhibiting the expression of E2F1, it is helpful to finally achieve the purpose of inhibiting the development of HCC.As our current understanding of E2F1 is insufficient, we sought to study the mechanism and related regulation of E2F1 in HCC. Thus, based on bioinformatics, we adopted the ceRNA mechanism as an entry point to analyze the lncRNAs and miRNAs related to E2F1 and explore the correlation between E2F1 and immune infiltration levels of various types of immune cells in HCC to identify potential biological targets and to prepare for subsequent basic research. We present the following article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-22-674/rc).

Methods

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). This research is divided into the following parts. First, we preliminarily verified that E2F1 is significantly overexpressed in malignant tumors including HCC through the TCGA database. Kaplan-Meier analysis was used to discuss the relationship between E2F1 expression and survival. Second, a retrospective analysis of 364 patients was performed to explore the relationship between E2F1 and clinicopathological parameters. The specificity and sensitivity of E2F1 as a prognostic indicator were also evaluated. Second, a retrospective analysis of 364 patients was performed to explore the relationship between E2F1 and clinicopathological parameters. The specificity and sensitivity of E2F1 as a prognostic indicator were also evaluated. Third, bioinformatics analysis was used to explore potential ceRNA mechanisms, and a total of 6 potential E2F1-related signaling pathways were screened. Then, the relationship between E2F1 and tumor-related immunity was explored using the Tumor IMmune Estimation Resource (TIMER) database. Finally, the drug sensitivity of E2F1 as a therapeutic target in HCC was explored.

Differential expression and survival analyses of E2F1

HCC-related transcriptome data and patient clinical information downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). A differential expression map of E2F1 in pan-cancer was obtained from the TIMER2.0 database (http://timer.cistrome.org/). Based on the magnitude of the P value of E2F1 for different cancer types, the P value was divided into P<0.001, P<0.01, and P<0.05. The differential expression analysis of E2F1 in HCC was analyzed and plotted using the “limma”, “ggplot2,” and “ggpubr” packages. Prognosis-related survival curves were downloaded from the Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn).

Analysis of clinical prognostic factors

The correlation between E2F1 expression and each prognostic factor was analyzed and plotted using the “limma” and “ggpubr” packages. Heatmaps for each clinical prognostic factor were drawn using the “limma” and “ComplexHeatmap” packages. The receiver operator characteristic (ROC) curves, calibration curves, and nomograms were made using the “survival,” “survminer,” “timeROC,” “regplot,” and “rms” analysis packages. Univariate and multivariate Cox analyses were conducted, and forest plots were generated using the “survival” package.

Establishment of mRNA-miRNA-lncRNA co-expression network and survival analysis of miRNA and lncRNA

The mRNA-miRNA and miRNA-LncRNA interaction data were downloaded from the starBase database (http://starbase.sysu.edu.cn/), with a programNum ≥2 as one of the mRNA-miRNA screening criteria. The correlation coefficient values (an R value >0.2 was defined as a positive correlation, and an R value <–0.2 was defined as negative correlation), differential expression values (a P value <0.01 was considered statistically significant), and survival curve values (a P value <0.05 was considered statistically significant) were screened out and plotted using the R language package. A conceptual diagram of the potential ceRNA mechanisms associated with E2F1 was drawn with BioRender.

E2F1 immune correlation analysis

The correlation analyses between E2F1 and various immune cells and immune checkpoints were visualized using various R language packages, including “limma,” “reshape2,” “ggplot2,” “ggpubr,” “vioplot,” “ggExtra,” and “corrplot.” The p values were calculated using the Spearman statistical method. A positive correlation was defined as a P value <0.05, an R value >0.2, a negative correlation was defined as a P value <0.05, an R value <–0.2, and a P value >0.05 was defined as not significant. E2F1, PDCD1, CD274, CTLA4, and LAG3 were analyzed using the TIMER 2.0 database (http://timer.cistrome.org/).

E2F1-related drug sensitivity evaluation

half maximal inhibitory concentration (IC50) represents the concentration required for the 50% inhibition of drug concentration. We calculated the IC50 of drugs using the “pRRophetic” R package with its dependencies “car, ridge preprocessCore, genefilter, and sva,” which contained information on the effects of 138 drugs. Boxplots were drawn using the “ggplot2” R package. A P value <0.05 indicated a statistically significant difference.

Statistical analysis

Wilcoxon rank-sum test was used to compare the difference between the two groups. Differential expression data were analyzed by “DESeq2” and “survival” R software. KM survival analysis was used for ROC curve analysis, univariate and multivariate Cox regression analysis. Spearman’s test was used to measure correlations between E2F1 and immune functions. And P value <0.05 was regarded as the significant threshold.

Results

Differential expression and survival analyses of E2F1 in HCC

We analyzed the expression of E2F1 in 38 cancer types in the TIMER2.0 database and found that E2F1 was significantly differentially expressed between the tumor group and the normal group in terms of 20 malignant tumors, including HCC (P<0.001; see ). We downloaded the HCC-related transcriptome data and patient clinical information downloaded from TCGA database, and found that E2F1 was significantly differentially expressed in malignant tumors (P<0.01; see ). We also searched E2F1-related disease-free survival (DFS) and overall survival (OS) of HCC patients in the GEPIA database and found that there were significant differences in the prognosis of the high-risk group (182 cases) and the low-risk group (182 cases) (DFS: P=0.0027, OS: P=0.0025), and the prognosis of the low-risk group was significantly better than that of the high-risk group in terms of both OS (see ) and DFS (see ).

Figure 1

Differential expression and prognostic survival curve of E2F1 (A) Differential expression of E2F1 in HCC from the TIMER2.0 database; (B) differential expression of E2F1 between the tumor group and the normal group in HCC; (C) pairwise differential analysis of E2F1 in HCC (*, P<0.05; **, P<0.01; ***, P<0.001); (D) overall survival curve from the GEPIA database; (E) disease-free survival curve from the GEPIA database. HCC, hepatocellular carcinoma; TIMER, Tumor IMmune Estimation Resource; GEPIA, Gene Expression Profiling Interactive Analysis.

Analysis of prognostic factors of E2F1 in HCC

Subsequently, we analyzed the key prognostic factors related to E2F1 and found that there was a significant difference between stage I, stage II, and stage III, and between stage III and IV (P<0.01; see ). Significant differences were found among all the grades, except grade 3 and grade 4 (see ). The expression of E2F1 differed significantly between T1 and T2, T3, and T4 (P<0.01; see ), and there were significant differences in tumor (T) stage, stage, and grade in the high- and low-expression groups (P<0.001; see ). In the sensitivity and specificity analyses, the ROC curve of the target gene E2F1 showed that the areas under the curve (AUCs) at 1, 3, and 5 years were 0.646, 0.628, and 0.584, respectively (see ). A nomogram was drawn to assess whether E2F1 could predict survival time in HCC (see ), and the feasibility of this prediction method was validated with a calibration curve (see ). Finally, we concluded that E2F1 expression and HCC stage were independent risk factors for prognosis through univariate and multivariate Cox regression analyses (see ).

Figure 2

E2F1 prognostic factors (A) stage; (B) grade; (C) T stage; (D) heat map; E2F1 was found to be correlated with diagnosis and prognosis; (E) calibration curve (OS); (F) nomogram; (G) the predictive effect of E2F1 (AUC); (H) univariate Cox regression analysis of prognosis-related risk factors; and (I) multivariate Cox regression analysis of prognosis-related risk factors. **, P<0.01; ***, P<0.001. OS, overall survival; AUC, area under curve.

Establishment of mRNA-miRNA co-expression network and related miRNA survival analysis

After downloading the E2F1-miRNA interaction data from the starBase database and using the R language package for the analysis, we screened 2 groups of miRNAs that were co-expressed and negatively correlated with E2F1 (i.e., had a correlation coefficient

Figure 3

Correlation, difference, and survival curve analyses between E2F1 and 2 miRNAs. (A) miR-29b-3p; (B) miR-29c-3p.

Next, we downloaded the lncRNA data that interacted with miR-29b-3p and miR-29c-3p from the starBase database, and screened and analyzed the correlation coefficients (those with a correlation coefficient value >0.2), log fold change (FC) values (those with a log FC value >0), survival curves, and the differential expression between the tumor group and the normal group (P<0.01) by R language. We also selected lncRNAs whose expression levels were positively correlated with E2F1 according to the above screening results (those with a correlation coefficient value >0.2, and a P value <0.01). Finally, the miR-29b-3p-related lncRNAs (i.e., LINC01224, PCBP1-AS1, and ITGA9-AS1), and the miR-29c-3p-related lncRNAs (i.e., SNHG7, THUMPD3-AS1, and LINC02323) were screened (see ). Based on the above results, the possible potential ceRNA mechanism diagram for E2F1 was constructed (see ).

Figure 4

Correlation, difference, and survival curve analyses between miR-29b-3p-E2F1-related lncRNA and miRNA and E2F1, respectively (A). (B) Correlation, difference, and survival curve analyses between miR-29c-3p-E2F1-related lncRNA and miRNA and E2F1, respectively.

Figure 5

Conceptual diagram of the E2F1-related ceRNA mechanism in HCC. HCC, hepatocellular carcinoma.

Correlation analyses of E2F1 with various immune cells and immune checkpoints

Additionally, we analyzed the correlations between immune cells and the levels of immune infiltration for E2F1, and found that cluster of differentiation (CD)4 memory activated T cells, memory B cells, eosinophils, and follicular helper T cells were positively correlated with E2F1 (R>0.2, P<0.01), and monocytes, naïve B cells, and CD4 memory resting T cells were negatively correlated with E2F1 (R<–0.2, P<0.01; see ). The TIMER database-related immune checkpoint analysis showed that E2F1 was positively correlated with PDCD1, CTLA4, and LAG3 (R>0.2, P<0.01; see ).

Figure 6

Correlations between E2F1 and various immune cells and immune checkpoints.

Drug sensitivity evaluation

We examined the relationship between the risk score and the IC50 of various drugs used in the clinical treatment of HCC, including imatinib, etoposide, and paclitaxel. Patients in the high-expression group appeared to be more susceptible to most drugs than those in the low-expression group (P<0.001; see ).

Figure 7

The IC50 of 39 drugs in E2F1 high- and low-expression groups (P<0.001).

Discussion

E2F1 was the first member of the E2F transcription factor family, which comprises 8 proteins, to be discovered (14). Based on their different functions, they are usually classified as activators (E2F1-e2f3a) or inhibitors (E2F3b-E2F8) (15). Studies have shown that E2F1 mainly regulates the transcription of S-phase cyclins and related genes required for DNA replication, DNA repair, and apoptosis (16). At present, the common genes that cause the abnormal activation of E2F1 mainly include retinoblastoma (Rb), Ras, and PI3K. The abnormal activation of E2F1 affects its downstream transcriptional targets, resulting in DNA replication stress (16). Its transcriptional targets include cyclin E and RRM2. Cyclin E promotes the phosphorylation of essential DNA replication factors to initiate and allow the progression of bidirectional DNA synthesis. Cyclin E overexpression results in enhanced CDK2 activity and cell cycle progression, thereby reducing the ability of cells to regulate the G1 (DNA prophase)-S (DNA replication period) transition (17). This regulatory mechanism has been widely observed in a number of malignancies (17-19). In addition, another important transcriptional target of E2F1 that could contribute to DNA replication stress is RRM2 (20), and the above signaling pathway of E2F1 has been reported in adrenocortical carcinoma (21), colorectal cancer (22), pancreas cancer (23), and other malignant tumors.

Using the TIMER database, we sought to identify the immune cells correlated with E2F1 in terms of the level of immune infiltration in HCC. We found that the expression of E2F1 was positively correlated with CD4 memory activated T cells, memory B cells, eosinophils, and follicular helper T cells, and negatively correlated with monocytes, naïve B cells, and CD4 memory resting T cells. Studies have shown that HTLV-1 basic leucine zipper factor (HBZ) is a related viral factor required for the viral replication and transformation of infected cells. HBZ protein interacts with the Rb/E2F-1 complex and induces the transcription of E2F target genes. The activation of the Rb/E2F pathway by the HBZ protein accelerates G1/S transition and apoptosis in primary CD4+ T cells (24). The downregulation of E2F1 decreases the susceptibility of CD8+ T cells. E2F1 has been shown to be a transcription factor for TBX21, a Th1 cell-specific transcription factor that controls the expression of the hallmark Th1 cytokine and interferon gamma (IFN-γ) (25). Thus, E2F1 plays an important role in tumor immunity by affecting the activation of effector CD8+ T cells (26).

E2F1 also significantly represses the transcriptional activity of the interleukin (IL)-6 promoter, while the overexpression of E2F1 promotes this activity. E2F1 regulates macrophage cytokine expression via IL-6 in nasopharyngeal carcinoma (NPC) cell supernatants, which supports its utility in the tumor microenvironment (TME). In a xenograft tumorigenesis model, small interfering–RNAs targeting E2F1 or E2F3 significantly inhibited tumor growth and reduced immune cell infiltration in the TME (27), which suggests that E2F1 can be regulated by modulating macrophage function. Further, E2F1 transactivates the IL-6 promoter, a very important inflammatory cytokine. However, E2F1 mostly acts as an inhibitor to negatively regulate dendritic cells (28), but its activation in mouse bone marrow-derived dendritic cells (DC2.4) cells is decreased by E2F1 knockdown and enhanced by E2F1 overexpression. The mechanism underlying this phenomenon is unclear; however, it may be related to the activation of p38 mitogen-activated protein kinase (MAPK) by E2F1, which directly promotes the activation of DC2.4 cells.

According to previous findings, the silencing of LINC01224 downregulates CHEK1 expression by competitively binding to miR-330-5p, thereby inhibiting HCC progression. Additionally, LINC01224 has been shown to induce HCC progression in vitro and accelerate HCC formation in nude mice by increasing CHEK1 expression (29). There are differences in the expression of PCBP1-AS1 in HCC. Notably, PCBP1-AS1 promotes HCC progression and HCC cell metastasis by combining with PCBP1 and regulating the PCBP1/PRL-3/serine/threonine kinase (AKT) pathway (30). The expression of lncRNA SNHG7 is upregulated in HCC, and elevated SNHG7 expression is closely associated to the staging, grading, vascular invasion, and poor prognosis in HCC patients. SNHG7 promotes HCC progression by regulating miR-122-5p and RPL4 (31). Additionally, studies have confirmed that low expression of miR-29b-3p, miR-29c-3p is associated with tumor growth, multiple pathological features, and shorter OS (32). Several HCC-related reports have noted that the overexpression of miR-29b-3p, miR-29c-3p significantly inhibits the proliferation, apoptosis, migration, and tumor growth of HCC cells in vivo (33,34).

In this study, using TCGA, GEPIA and starBase databases, we identified the miRNAs (i.e., miR-29b-3p and miR-29c-3p) related to transcription factor E2F1 in HCC by R language. We also used the lncRNAs related to E2F1 (i.e., miR-29b-3p, and miR-29c-3p) to construct ceRNA models. Further, we analyzed the related immune cell infiltration, immune checkpoints, and drug sensitivity of E2F1 using the TIMER database. It should be noted that this research was based on a bioinformatics analysis; thus, the validity of the findings needs to be further verified by basic experimental research. However, our results still provide a very valuable direction and reference for research on transcription factor E2F1, which may be helpful in identifying research targets for future HCC-related molecular biological therapy and immunotherapy.

The article’s supplementary files as