The data on the pan-carcinoma project were downloaded from TCGA database. With the use of the R language (version 3.63), the Mann–Whitney U test was utilized to examine the difference in expression among cancer and normal tissues, and the ggplot2 package was utilized for visualization.

In this study, immunohistochemical images of E2F1, E2F2, and E2F8 in THCA and normal tissues were identified and extracted from HPA (https://www.proteinatlas.org/) for comparison (Thul et al., 2017).

We used the GEPIA2 (http://gepia2.cancer-pku.cn/#analysis) database to study the expression profile of E2Fs in THCA and normal thyroid tissues and to investigate the association between their level and clinical stages of THCA patients (Tang et al., 2019).

The UALCAN database (http://ualcan.path.uab.edu/) is an established secondary site based on the analysis of TCGA database (Chandrashekar et al., 2017). In this study, the gene name of E2F family was first input and the cancer type was selected as thyroid cancer to study the level of the E2F gene family in THCA and its connection with clinicopathological features of THCA patients.

Appropriate THCA datasets were selected from the cBioPortal database (https://www.cbioportal.org/) to analyze the mutation rate of the E2F gene family in THCA (Cerami et al., 2012), to explore whether the expression changes of the E2F gene in THCA are caused by gene mutations, and to speculate the expression stability of E2Fs in THCA.

The TCGA thyroid cancer project’s Level 3 HTSeq-FPKM format RNAseq data were analyzed using the R language (version 3.63) to examine the expression association between eight E2F family members using Spearman’s correlation analysis, and the ggplot2 package was used for visualization.

The STRING database (http://stringdb.org) is used to assess and generalize data based on protein–protein interactions, including their functional and physical correlations (Szklarczyk et al., 2017). The E2F gene was input into the homepage of the official website of the GeneMANIA database (https://genemania.org/) to construct a dendrogram of its related gene network (Franz et al., 2018).

Using the online analytical tool VENN (http://bioinformatics.psb.ugent.be/webtools/Venn/), Venn diagrams were constructed to obtain the cross genes shared by the two gene datasets (Bardou et al., 2014).

In this research, we used the GSCA (http://bioinfo.life.hust.edu.cn/web/GSCALite/) online tool to investigate the association between prognosis and E2F family in various tumors (Liu et al., 2018).

The patients’ clinical information of 510 THCA tissues and 58 normal tissues were obtained from TCGA database. The clinical information of patients with THCA was integrated, and the Kaplan–Meier analysis with the log-rank test was performed using the survival package (version 3.2-10). The results from the survival analysis were visualized using the survminer package (version 0.4.9).

In this study, the LinkInterpreter module of LinkedOmics (http://www.linkedomics.org/login.php) was used to perform ORA to explore the signaling pathways associated with the participation of the E2F family in THCA. Pathways with FDR ≤ 0.05 were considered significant (Vasaikar et al., 2018).

According to the median gene expression, tumor samples from TCGA project were divided into high expression and low expression groups. To identify core gene-related signaling pathways, c2.cp.kegg.v7.2.symbols.gmt was chosen as a reference using GSEA software for analysis (Subramanian et al., 2005). By examining 1000 permutations, the normalized enrichment analysis score (NES) was calculated. According to the previous research, |NES| ≥ 2.0, p < 0.05, and false discovery rate (FDR) < 0.25 were used as thresholds according to the standard screening gene set (Fu et al., 2022).

The TISIDB database (http://cis.hku.hk/TISIDB/) is one of the commonly used databases for cancer immune infiltration markers (Ru et al., 2019). E2F family members were input into the TISIDB database, and the association between their expression and that of immune checkpoint inhibitor genes was examined.

The TIMER database (https://cistrome.shinyapps.io/timer/) is one of the commonly used cancer immune infiltration marker databases (Li et al., 2017). The correlation between E2F family and immune cells was explored using the gene module in the TIMER database. Subsequently, we used the somatic copy number alteration (SCNA) module to explore the relationship between tumor infiltration level and somatic copy number changes of the E2F family in THCA. Finally, the correlation analysis module was employed to investigate the association between E2F levels and immune-infiltrating cell markers.

Survminer, survival, and ggplot2 packages of R software (version 3.6.3) were utilized. The Mann–Whitney U test was utilized for differential expression analysis. Kaplan–Meier and log-rank sum tests were utilized for prognostic analysis. The statistical significance was set at p < 0.05.

To study the levels of the E2F family in pan-carcinoma, RNAseq data from TCGA were obtained, and the results were analyzed and calculated by the R language and visualized by the ggplot2 package. The outcomes displayed that the E2F gene was greatly expressed in most cancers (Figures 1, 2). For THCA, E2F1, 2, 3, 7, and 8 genes were strongly expressed, indicating that the E2F family plays a significant carcinogenic role in THCA (Figures 1, 2). Next, we investigated the protein level of E2Fs in THCA. Immunohistochemical results of three E2F members, namely, E2F1, E2F2, and E2F8, in THCA were found in the HPA database. The outcomes displayed that the protein level of E2F1 was higher in THCA, while E2F2 and E2F8 did not significantly vary from normal tissues in terms of protein expression levels (Supplementary Figure S1). The findings of the analysis of the association between E2F transcription levels and THCA tumor stages using the GEPIA2 database displayed that the level of E2F2 and E2F5 was associated with the THCA cancer stage (Supplementary Figure S2). Finally, we discovered that the E2F gene family was related to race, gender, age, TNM stage, pathological stage, histological type, extrathyroidal extension, and primary neoplasm focus type in THCA (Supplementary Table S1).

The genetic alterations of E2Fs in the TCGA-THCA samples were investigated using the cBioPortal database. The findings demonstrated that the E2F family was highly conserved in THCA and rarely mutated (Figure 3A). We then used the muTarget database to explore the expression of genes that is affected by the E2F mutation in thyroid cancer and genes themselves having mutations that affected the expression of E2Fs. Finally, only mutations in NRAS in thyroid cancer were found to reduce E2F8 expression. This is consistent with the highly conserved results of E2Fs in thyroid cancer. Then, we used the R language to perform correlation analysis to explore the expression association among the eight members of the E2F family in THCA. The ggplot2 package was used for visualization. The findings demonstrated that the vast majority of E2F family members were positively correlated (Figure 3B). TP53, a key gene of the apoptosis pathway, and AKT1, a key gene of the cell survival pathway, interact closely with E2Fs. Therefore, we explored the expression patterns between these markers and E2Fs. The results of the chord diagram showed that there was a significant positive correlation between them (Supplementary Figure S3). Finally, the gene and protein networks that potentially interacted with the E2F family were investigated using the GeneMANIA and STRING databases. Ten target proteins and 20 target genes were identified to interact with the E2F family (Figures 3C, D). Cross-analysis of the aforementioned two groups yielded a common molecule, TFDP1 (Figure 3E).

We divided tumor cases into two groups based on the level of expression of each E2F members; next, we utilized TCGA dataset to study the relationship between E2F gene level and THCA patients’ prognoses. A pan-cancer study utilizing the GSCA online tool revealed, as shown in Figure 4A, that the E2F family was associated with the prognosis of a range of malignancies and that the stronger the E2F family expression, the worst the prognosis. Next, the survival package was used to analyze the prognosis of RNAseq data of the TCGA-THCA project. As THCA has a better prognosis, a progression-free interval (PFI) was chosen instead of overall survival (OS) as the outcome event. The findings demonstrated that patients with high E2F2, E2F7, and E2F8 expression had shorter PFI (Figure 4B). The subgroup K–M analysis was then used to assess the association between the level of E2F members and clinicopathology of THCA in order to further understand the association between E2F family and clinical characteristics of THCA. We discovered that the high E2F1 level was associated with poor PFI in patients with THCA in age, T stage, and pathological stage (Supplementary Table S2). Finally, to investigate the independent risk variables of PFI in THCA, regression analyses were carried out. Univariate Cox analysis showed that T stage, M stage, pathologic stage, extrathyroidal extension, E2F2, E2F7, and E2F8 were the risk factors for PFI (p = 0.001, 0.001, 0.001, 0.023, 0.017, 0.025, and 0.046, respectively). Only M stage (p = 0.032), according to multivariate Cox regression analysis, was identified as an independent risk factor for PFI (Supplementary Table S3).

We conducted ORA utilizing the LinkedOmics database to examine the function of E2Fs in THCA. The outcomes displayed that the F2F family was closely associated with Ras signaling pathway, PDGF signaling pathway, apoptosis, and immune response (Figure 5). The signaling pathways impacted by the expression of the E2F family in THCA were then investigated using single-gene GSEA. The findings showed that the functions of the E2F family may focus on GPCR ligand binding, neutrophil degranulation, and signaling by interleukins (Figure 6).

By regulating the activity of related immune cells, the immune checkpoint inhibitor (ICI) improves anti-tumor immune response, which is a major progress in tumor treatment in recent years. Therefore, the TISIDB database was used to explore whether E2Fs in THCA is associated with immune checkpoint inhibitors, for the purpose of exploring whether the function of the E2F family in THCA development is related to the immunosuppressive effect. We included all eight members of the E2F family in this study, and the findings demonstrated that their expression levels correlated differently with each immune checkpoint inhibitor (Figures 7A–H).

Immune infiltration and cancer microenvironment develop a significant function in the emergence and development of various malignancies. Different clinical outcomes might result from the dynamic interplay of tumor genetic features and tumor microenvironment. We studied the association between the E2F family and tumor-infiltrating lymphocytes (TILs) in THCA. The findings demonstrated that the E2F family members were strongly related to TIL, among which B cells and CD8⁺ T cells were positively correlated with all E2F family members. CD4⁺ T cells, macrophages, and neutrophils were positively correlated with E2F family members except E2F1. Dendritic cells were positively correlated with E2F family members except E2F1 and E2F6 (Figures 8A–H). The SCNA template analysis in the TIMER database showed that the E2F family members’ arm-level deletion and arm-level gain in THCA were strongly related to the level of TIL (Figure 9). Then, we assessed the association between E2F level and numerous immune infiltration-related indicators in order to further investigate the link between tumor-infiltrating immune cells and E2F family members in THCA. The results showed a great association between the level of E2F2, 7, and 8 and most of the ensembles of immune-infiltrating cell markers (Supplementary Table S4). In particular, for the T-cell exhaustion-associated marker ensemble, the outcomes were in line with Figure 7.