# Analysis of Transcription Factor-Related Regulatory Networks Based on Bioinformatics Analysis and Validation in Hepatocellular Carcinoma

Authors: Shui Liu; Xiaoxiao Yao; Dan Zhang; Jiyao Sheng; Xin Wen; Qingyu Wang; Gaoyang Chen; Zhaoyan Li; Zhenwu Du; Xuewen Zhang
Journal: BioMed Research International
 (2018)

---

## Abstract

Hepatocellular carcinoma (HCC) accounts for a significant proportion of liver cancer, which has become the second most common cause of cancer-related mortality worldwide. To investigate the potential mechanisms of invasion and progression of HCC, bioinformatics analysis and validation by qRT-PCR were performed. We found 237 differentially expressed genes (DEGs) including EGR1, FOS, and FOSB, which were three cancer-related transcription factors. Subsequently, we constructed TF-gene network and miRNA-TF-mRNA network based on data obtained from mRNA and miRNA expression profiles for analysis of HCC. We found that 42 key genes from the TF-gene network including EGR1, FOS, and FOSB were most enriched in the p53 signaling pathway. The qRT-PCR data confirmed that mRNA levels of EGR1, FOS, and FOSB all were decreased in HCC tissues. In addition, we confirmed that the mRNA levels of CCNB1, CCNB2, and CHEK1, three key markers of the p53 signaling pathway, were all increased in HCC tissues by bioinformatics analysis and qRT-PCR validation. Therefore, we speculated that miR-181a-5p, which was upregulated in HCC tissues, could regulate FOS and EGR1 to promote the invasion and progression of HCC by p53 signaling pathway. Overall, the study provides support for the possible mechanisms of progression in HCC.


---

## Body

## 1. Introduction

Hepatocellular carcinoma (HCC) is one of the common digestive system malignancies with high mortality, which has become the second most common cause of cancer-related mortality worldwide [1, 2]. There are more than 466,100 estimated new cases and 422,100 estimated death cases every year in China [3]. Although there is extensive research about the therapies for HCC, the specific mechanisms of HCC occurrence and development were not clear. It is significant to investigate the underlying mechanisms of HCC invasion and progression to develop effective diagnostic and therapeutic strategy.Over the last decades, bioinformatics and microarray technology have been widely used to screen the molecular mechanisms of HCC, which provide powerful research support for the diagnosis and treatment of HCC [4–6]. By bioinformatics analysis and microarray technology, functions of some differentially expressed genes (DEGs) in HCC have been explored, which paved the way for exploring complicated process in the occurrence and development of HCC [5, 7–9]. As an important part of participating in life activities, transcription factors (TFs) have been reported to play an important role in the development of HCC by numerous studies [10, 11]. A work by S. Hua et al. found that ETS1, a cancer-related TF, could through interaction with miR-139-5p inhibit cell proliferation, migration, and invasion in HCC [12]. Another research reported that NFκB/EGR1 signaling pathway induced miR-3928v to promote the progression of HCC [13]. Even so, more research is still needed to explore the specific role of different TFs in the progression of HCC.During the present study, we analyzed 4 mRNA microarray datasets and 1 miRNA microarray dataset from Gene Expression Omnibus (GEO) to obtain DEGs and differentially expressed miRNAs (DEMs) between HCC tissues and adjacent noncancerous tissues by GEO2R. Subsequently, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, and protein-protein interaction (PPI) network analysis were performed to reveal the interaction relationship between DEGs and DEMs and to explore the underlying molecular mechanisms in the carcinogenesis and progression of HCC [14–17]. In addition, we screened differentially expressed TFs in HCC, constructed miRNA-TF-mRNA networks, and proposed a potential miRNA-TF-signaling pathway axis, which identified a systematic exposition of the relevant transcriptional regulation modes associated with invasion and progression of HCC. This study could provide insightful information and the valuable clue for biomarker discovery and novel treatment strategy in HCC.


## 2. Materials and Methods

### 2.1. Microarray Data Information and DEGs Identification

NCBI-GEO is a free database of microarray/gene profile and next-generation sequencing, from which HCC and adjacent nontumor tissue gene expression profile of GSE84402, GSE33006, GSE84005, GSE12941, and GSE64632 were obtained [4, 18, 19, 21]. The DEGs and DEMs between HCC tissues and adjacent nontumor tissues were screened using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r), which was an interactive web tool for identifying DEGs across experimental conditions among two or more datasets in a GEO series. Identification of commonly differentially expressed mRNAs from the four cohort profile data sets (GSE84402, GSE33006, GSE84005, and GSE12941) was performed by FunRich software (available online: http://www.funrich.org/). |logFC| (fold change) >1 and P value <0.05 were considered statistically significant. The pipeline of the whole process of this study was shown in Figure 1.Figure 1
Process of TF-related regulatory network construction and key genes identification in HCC.


### 2.2. Construction of Transcription Factor Networks

Transcriptional regulatory element database (TRED) (https://cb.utdallas.edu/cgibin/TRED/tred.cgi?process=home) is a database of transcriptional regulatory elements, which provides 36 cancer-related TFs and corresponding regulatory networks [22]. Transcriptional Regulatory Relationships Unraveled by Sentenced-based Text mining version 2.0 database (TRRUST) (available online: http://www.grnpedia.org/trrust/) is a database based on the literature to reflect the relationship between transcriptional regulation [23]. And the prediction of miRNA-target genes was performed by miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) [24]. According to the regulatory interaction, TF-target network and miRNA-TF-mRNA network were constructed based on gene expression profile, TRRUST version 2, TRED, and miRTarBase by Cytoscape 3.6.0 software.


### 2.3. Functional and Pathway Enrichment Analysis

Omicsbean is an online biological information database that integrates biological data and analysis tools and provides a comprehensive set of functional annotation information of genes and proteins for users to extract biological information. To analyze the function of DEGs, GO and KEGG pathway enrichment analysis were performed using Omicsbean online database (http://www.omicsbean.com:88). P value<0.05 was considered statistically significant.


### 2.4. Validation in TCGA Dataset and Modular Analysis of the Key Genes

Validation of the key genes in the TCGA Dataset was performed by UALCAN (http://ualcan.path.uab.edu/index.html) [25]. We performed Kaplan-Meier plots and boxplots to verify the differential expression of the key genes between primary tumor and normal liver and effect of key genes expression level on LIHC patients' survival. Correlation analysis was performed by Linkedomics (http://www.linkedomics.org) [26].


### 2.5. Validation of the Key Genes by Quantitative Real-Time RT-PCR (qRT-PCR)

A total of 20 HCC patients were recruited for tumor and adjacent nontumor tissues collection from the China-Japan Union Hospital of Jilin University, Changchun, China. This study was approved by the Ethics Committee of the Second Clinical Medical College, Jilin University, and each patient consented in a written informed consent form. The data were analyzed anonymously. All tissues were taken from the surgery room and snap-frozen and stored in liquid nitrogen within 10 min after the resection. The clinicopathological characteristics of 20 HCC patients were shown in TableS1.Tissue RNA was isolated using Trizol (Invitrogen, CA, USA) and further purified using the MiniBEST Universal RNA Extraction Kit (TaKaRa, China) according to the manufacturer's instructions. RNA concentration was then measured using the NanoDrop 2000 spectrophotometer (Thermo Scientific, MA, USA) with A260/A280 ratio in the range 1.8~2.0 and RNA concentration ranged from 100 ng/µl to 1 µg/µl.For qRT-PCR analysis, less than 5µg total RNA including microRNA was reverse transcribed to cDNA with 1st strand cDNA Synthesis Kit (Takara, China) and miRNA First Strand cDNA Synthesis (Sangon, China); the expression of mRNA and microRNA was examined by qRT-PCR with TransStart ® Top Green qPCR SuperMix (TransGen Biotech, China), microRNAs qPCR Kit (Sangon, China), and Applied Biosystems 7500 Fast Real-Time PCR System. The relative expression of mRNA and microRNA was normalized to GAPDH/U6 expression by comparative Ct method (2-ΔΔCt, ΔCt =Cttarget-CtGAPDH/U6, and ΔΔCt=ΔCttumor-ΔCtnon-tumor). All primers were designed with Primer 7.0 Software; primer sequences for amplification were listed in Table 8.


### 2.6. Statistical Analysis

Data from qRT-PCR were analyzed with GraphPad Prism version 7.0, and differences between the two groups were statistically evaluated by two-tailed Student’s t-test with p value <0.05 considered as significant. The Pearson correlation coefficient was used to examine the relationship between key genes.


## 2.1. Microarray Data Information and DEGs Identification

NCBI-GEO is a free database of microarray/gene profile and next-generation sequencing, from which HCC and adjacent nontumor tissue gene expression profile of GSE84402, GSE33006, GSE84005, GSE12941, and GSE64632 were obtained [4, 18, 19, 21]. The DEGs and DEMs between HCC tissues and adjacent nontumor tissues were screened using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r), which was an interactive web tool for identifying DEGs across experimental conditions among two or more datasets in a GEO series. Identification of commonly differentially expressed mRNAs from the four cohort profile data sets (GSE84402, GSE33006, GSE84005, and GSE12941) was performed by FunRich software (available online: http://www.funrich.org/). |logFC| (fold change) >1 and P value <0.05 were considered statistically significant. The pipeline of the whole process of this study was shown in Figure 1.Figure 1
Process of TF-related regulatory network construction and key genes identification in HCC.


## 2.2. Construction of Transcription Factor Networks

Transcriptional regulatory element database (TRED) (https://cb.utdallas.edu/cgibin/TRED/tred.cgi?process=home) is a database of transcriptional regulatory elements, which provides 36 cancer-related TFs and corresponding regulatory networks [22]. Transcriptional Regulatory Relationships Unraveled by Sentenced-based Text mining version 2.0 database (TRRUST) (available online: http://www.grnpedia.org/trrust/) is a database based on the literature to reflect the relationship between transcriptional regulation [23]. And the prediction of miRNA-target genes was performed by miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) [24]. According to the regulatory interaction, TF-target network and miRNA-TF-mRNA network were constructed based on gene expression profile, TRRUST version 2, TRED, and miRTarBase by Cytoscape 3.6.0 software.


## 2.3. Functional and Pathway Enrichment Analysis

Omicsbean is an online biological information database that integrates biological data and analysis tools and provides a comprehensive set of functional annotation information of genes and proteins for users to extract biological information. To analyze the function of DEGs, GO and KEGG pathway enrichment analysis were performed using Omicsbean online database (http://www.omicsbean.com:88). P value<0.05 was considered statistically significant.


## 2.4. Validation in TCGA Dataset and Modular Analysis of the Key Genes

Validation of the key genes in the TCGA Dataset was performed by UALCAN (http://ualcan.path.uab.edu/index.html) [25]. We performed Kaplan-Meier plots and boxplots to verify the differential expression of the key genes between primary tumor and normal liver and effect of key genes expression level on LIHC patients' survival. Correlation analysis was performed by Linkedomics (http://www.linkedomics.org) [26].


## 2.5. Validation of the Key Genes by Quantitative Real-Time RT-PCR (qRT-PCR)

A total of 20 HCC patients were recruited for tumor and adjacent nontumor tissues collection from the China-Japan Union Hospital of Jilin University, Changchun, China. This study was approved by the Ethics Committee of the Second Clinical Medical College, Jilin University, and each patient consented in a written informed consent form. The data were analyzed anonymously. All tissues were taken from the surgery room and snap-frozen and stored in liquid nitrogen within 10 min after the resection. The clinicopathological characteristics of 20 HCC patients were shown in TableS1.Tissue RNA was isolated using Trizol (Invitrogen, CA, USA) and further purified using the MiniBEST Universal RNA Extraction Kit (TaKaRa, China) according to the manufacturer's instructions. RNA concentration was then measured using the NanoDrop 2000 spectrophotometer (Thermo Scientific, MA, USA) with A260/A280 ratio in the range 1.8~2.0 and RNA concentration ranged from 100 ng/µl to 1 µg/µl.For qRT-PCR analysis, less than 5µg total RNA including microRNA was reverse transcribed to cDNA with 1st strand cDNA Synthesis Kit (Takara, China) and miRNA First Strand cDNA Synthesis (Sangon, China); the expression of mRNA and microRNA was examined by qRT-PCR with TransStart ® Top Green qPCR SuperMix (TransGen Biotech, China), microRNAs qPCR Kit (Sangon, China), and Applied Biosystems 7500 Fast Real-Time PCR System. The relative expression of mRNA and microRNA was normalized to GAPDH/U6 expression by comparative Ct method (2-ΔΔCt, ΔCt =Cttarget-CtGAPDH/U6, and ΔΔCt=ΔCttumor-ΔCtnon-tumor). All primers were designed with Primer 7.0 Software; primer sequences for amplification were listed in Table 8.


## 2.6. Statistical Analysis

Data from qRT-PCR were analyzed with GraphPad Prism version 7.0, and differences between the two groups were statistically evaluated by two-tailed Student’s t-test with p value <0.05 considered as significant. The Pearson correlation coefficient was used to examine the relationship between key genes.


## 3. Results

### 3.1. Identification and Enrichment Analysis of DEGs in HCC

HCC and adjacent nontumor tissue gene expression profile of GSE84402, GSE33006, GSE84005, and GSE12941 were obtained from NCBI-GEO. We used p<0.05 and |logFC|>1 as cut-off criterion. After integrated bioinformatical analysis, a total of 237 differentially expressed genes were identified from the four profile datasets, including 57 upregulated genes and 180 downregulated genes (Tables1–3, Figure 2).Table 1
The information of expression profiles.
Series Platform  Tissues Adjacent tissues Tumor tissues Reference GSE84402 GPL570 mRNA Liver 14 14 [4] GSE33006 GPL570 mRNA Liver 3 3 [18] GSE84005 GPL5175 mRNA Liver 38 38 - GSE12941 GPL5175 mRNA Liver 10 10 [19] GSE64632 GPL18116 miRNA Liver 3 3 [20]Table 2
The differentially expressed genes of downloaded expression profiles (P value<0.05, FC>2).
Upregulation Downregulation Total GSE84402 2232 1363 3595 GSE33006 1680 2696 4376 GSE84005 481 530 1011 GSE12941 175 385 560 Co-overexpression DEGs 57 180 237Table 3
The common DEGs of four gene expression profiles (P value<0.05, FC>2).
Common DEGs Gene symbol Downregulated DEGs CLEC4M, CXCL14, NAT2, FOSB, CFP, CRHBP, MARCO, CYP1A2, CLEC4G, FCN3, CNDP1, ABCA8, APOF, GYS2, GLS2, BMPER, ADRA1A, OIT3, VIPR1, GBA3, CDA, SDS, MT1G, SHBG, CYP3A4, HAMP, HAO2, GREM2, CYP39A1, HGFAC, C7, ANGPTL1, TIMD4, DPT, FOS, INMT, NNMT, TAT, C9, DNASE1L3, AKR1D1, IGJ, CYP8B1, PGLYRP2, TDO2, GPM6A, HHIP, BBOX1, CYP2B6, MT1X, RDH16, NPY1R, ZG16, DCN, GNMT, F9, C8A, HPGD, CETP, MFAP4, PLAC8, SRD5A2, AFM, ITIH4, LIFR, STAB2, HGF, C14ORF68, SRPX, TMEM27, ASPA, KLKB1, EGR1, ALDH8A1, SAA4, PHYHD1, CYP2C9, C6, DIRAS3, GHR, GNA14, ALPL, CD5L, PCK1, STEAP4, BMP5, IGF1, AGXT2, PZP, PON3, ACMSD, EPHA3, HPX, PRG4, RBMS3, ST3GAL6, FBP1, HBB, HAL, FTCD, ENO3, COLEC10, STEAP3, PROZ, CXCL12, SPP2, PDE7B, P2RY13, SLC38A4, ECM1, C8ORF4, SLC7A2, PRELP, C8B, CCL19, PTPRB, ITGA9, CHST4,ENG, RAMP3, GDF2, GPR128, ABCA9, CCDC3, VSIG4, GZMK, RND3, MBL2, LAMA2, EMR1, SOCS2, F8, CCBE1, CD69, CD163, PTGIS, FXYD1, DNAJC12, CA2, UROC1, GCH1, WDR72, RASGEF1B, COLEC11, IGFBP3, RGN, IFITM1, FBLN5, CYP2C18, CD4, GCDH, PKHD1, OLFML3, FGFR2, TEK, CPS1, TRIM22, ID1, GNAO1, SLC5A1, DHODH, PLSCR4, THBS1, ADH6, VNN1, HSD17B2, AXL, ZFP36, BGN, PDCD1LG2, SLC17A2, ASPN, SERPINA4, EHD3, FXYD6, IGFBP1, ACADSB, ADAMTS1, SFRP5, C1ORF162 Upregulated DEGs ENAH, TXNRD1, MCM4, GNPAT, FAM38B, SORT1, SRXN1, ECT2, TARBP1, RGS5, RFX5, TKT, RAP2A, TP53I3, TUBG1, ITGA6, PRKDC, AKR1C3, CCT3, TLCD1, ATAD2, HSPB1, CENPF, NUSAP1, SLC44A5, SPP1, CHEK1, CDH13, PEG10, MKI67, ASPM, ROBO1, CCNA2, BUB1, CKAP2, CDC6, UCK2, DTL, AURKA, TPX2, MDK, DEPDC1, GPC3, CAP2, UBE2T, PTTG1, CDKN3, RRM2, PRC1, FOXM1, CCNB2, SPINK1, ANLN, KIF4A, PBK, CCNB1, GINS1Figure 2
Identification of 237 common DEGs from the four cohort profile data sets. (a) 57 commonly upregulated DEGs. (b) 180 commonly downregulated DEGs. Different color areas represented different datasets. The cross areas meant the commonly changed DEGs. DEGs were identified with classical t-test; statistically significant DEGs were defined with p<0.05 and |logFC|>1 as the cut-off criterion.
(a) (b)The enrichment analysis of DEGs was performed. As shown in Table4 and Figure 3, in the biological process group, the DEGs were mainly enriched in response to chemical. In the cellular component group, the DEGs were most enriched in the extracellular region. In the molecular function group, the DEGs were most enriched in protein binding. And the most enriched KEGG pathway was complement and coagulation cascades.Table 4
Enrichment analysis of DEGs.
(a)
Top 10 enriched biological processes of DEGs
GO_Name GO_ID P value Count response to chemical GO:0042221 8.60E-31 103 cellular response to chemical stimulus GO:0070887 9.67E-27 80 response to stimulus GO:0050896 2.83E-25 158 single-organism metabolic process GO:0044710 1.07E-24 129 response to organic substance GO:0010033 1.49E-22 75 small molecule metabolic process GO:0044281 2.85E-22 76 single-multicellular organism process GO:0044707 6.06E-20 109 negative regulation of biological process GO:0048519 6.44E-20 91 negative regulation of cellular process GO:0048523 8.35E-20 86 regulation of response to stimulus GO:0048583 9.64E-20 82  (b)
Top 10 enriched cell components of DEGs
GO_Name GO_ID P value Count extracellular region GO:0005576 1.17E-37 110 extracellular region part GO:0044421 3.22E-35 96 membrane-bounded vesicle GO:0031988 3.92E-29 83 extracellular space GO:0005615 5.47E-29 60 vesicle GO:0031982 5.72E-29 89 extracellular vesicle GO:1903561 7.65E-27 67 extracellular organelle GO:0043230 7.96E-27 67 extracellular membrane-bounded organelle GO:0065010 2.19E-26 66 extracellular exosome GO:0070062 3.55E-26 66 cytoplasmic part GO:0044444 6.80E-14 136  (c)
Top 10 enriched molecular functions of DEGs
GO_Name GO_ID P value Count protein binding GO:0005515 4.15E-21 146 identical protein binding GO:0042802 2.61E-09 30 protein homodimerization activity GO:0042803 3.74E-09 22 oxygen binding GO:0019825 1.97E-08 8 growth factor binding GO:0019838 4.52E-08 10 extracellular matrix binding GO:0050840 6.88E-08 7 oxidoreductase activity GO:0016491 7.54E-08 32 binding GO:0005488 1.17E-07 186 protein dimerization activity GO:0046983 1.97E-07 27 oxidoreductase activity, acting on CH or CH2 groups GO:0016725 2.42E-07 4  (d)
Top 10 enriched KEGG pathways of DEGs
Pathway Name Pathway ID P value Count Complement and coagulation cascades hsa04610 1.79E-06 10 Metabolic pathways hsa01100 1.78E-05 43 p53 signaling pathway hsa04115 3.88E-05 8 Prion diseases hsa05020 3.97E-05 6 Histidine metabolism hsa00340 6.78E-05 5 Retinol metabolism hsa00830 1.93E-04 7 Steroid hormone biosynthesis hsa00140 6.97E-04 6 Tryptophan metabolism hsa00380 8.26E-04 5 Cell cycle hsa04110 2.17E-03 8 Caffeine metabolism hsa00232 3.35E-03 2Figure 3
Enrichment analysis of DEGs.


### 3.2. Identification and Functional Analysis of Key Differentially Expressed Genes, Construction of TF-Related Networks in HCC

We intersected 237 common DEGs and curated 36 cancer-related TF families to get three differentially expressed TFs (EGR1, FOS, and FOSB) in HCC (Table5). We searched the target genes and regulators of the three TFs through TRRUST version 2 database. Based on the results, we found that STAT6, a cancer-related TF, regulated FOS and FOSB, while FOS and EGR1 were coregulated by 2 cancer-related TFs, including BRCA1 and SP1, respectively. Similarly, TP53 and NKFB1 were coregulated by FOS and EGR1, respectively. Since FOS and EGR1 are linked to each other by 4 cancer-related TFs including BRCA1, SP1, TP53, and NKFB1, we also included these four TFs in the subsequent data of network construction. We screened the 4 TFs in the TRRUST version 2 database query to get each TF target genes and the DEGs from four profile datasets in GEO and TF-target genes for intersection analysis to obtain 42 differentially expressed target genes of the TFs (key genes) (Table 6), which laid the foundation for our next step to construct the gene signaling regulatory network in HCC.Table 5
Curated 36 cancer-related TF families.
Family Full Name Members (Official Gene Symbols) AP1 Activator Protein 1 FOS, FOSB, JUN, JUNB, JUND  AP2 Activator Protein 2 TFAP2A, TFAP2B, TFAP2C, TFAP2D, TFAP2E  AR Androgen Receptor AR  ATF Activating Transcription Factor ATF1 – 7  BCL B-cell CLL/lymphoma BCL3, BCL6  BRCA breast cancer susceptibility protein BRCA1 – 3  CEBP CCAAT/enhancer binding protein CEBPA, CEBPB, CEBPD, CEBPE, CEBPG  CREB cAMP responsive element binding protein CREB1 - 5, CREM  E2F E2F transcription factor E2F1 – 7  EGR early growth response protein EGR1 – 4  ELK member of ETS oncogene family ELK1, ELK3, ELK4  ER Estrogen Receptor ESR1, ESR2  ERG ets-related gene ERG  ETS ETS-domain transcription factor ETS1, ETS2, ETV4, SPI1  FLI1 friend leukemia integration site1 FLI1  GLI glioma-associated oncogene homolog GLI1 - 4  HIF Hypoxia-inducible factor HIF1A, ARNT, EPAS1, HIF3A  HLF hepatic leukemia factor HLF  HOX homeobox gene HOXA, HOXB, HOXD series, CHX10, MSX1, MSX2, TLX1, PBX2  LEF lymphoid enhancing factor LEF1  MYB myeloblastosis oncogene MYB, MYBL1, MYBL2  MYC myelocytomatosis viral oncogene homolog MYC  NFI nuclear factor I; CCAAT-binding transcription factor NFIA, NFIB, NFIC, NFIX  NFKB Nuclear factor kappa B, reticuloendotheliosis oncogene NFKB1, NFKB2, RELA, RELB, REL  OCT Octamer binding proteins POU2F1 - 3, POU3F1 - 2, POU5F1  p53 P53 family TP53, TP73L, TP73  PAX paired box gene PAX1 – 9  PPAR Peroxisome proliferator-activated receptor PPARA, PPARD, PPARG  PR Progesterone Receptor PGR  RAR retinoic acid receptor RARA, RARB, RARG  SMAD Mothers Against Decapentaplegic homolog SMAD1 - 9  SP sequence-specific transcription factor SP1 – 8  STAT signal transducer and activator of transcription STAT1 - 6  TAL1 T-cell acute lymphocytic leukemia-1 protein TAL1  USF upstream stimulatory factor USF1, USF2  WT1 Wilms tumor 1 (zinc finger protein) WT1Table 6
Differentially expressed target genes of the transcription factor.
Gene symbol differentially expressed target genes of the transcription factor CXCL14, FOSB, CYP1A2, ADRA1A, CYP3A4, FOS, PGLYRP2, CYP2B6, CETP, HGF, EGR1, PCK1, IGF1, FBP1, CXCL12, CCL19, CHST4, ENG, F8, CD69, CD163, IGFBP3, RGN, TRIM22, ID1, THBS1, ADH6, HSD17B2, BGN, ASPN, IGFBP1, CCNB1, TP53I3, PRC1, MKI67, ASPM, FOXM1, PTTG1, CHEK1, CCNB2, CCNA2, SPP1Based on the above TRRUST version 2 database, we constructed the TF-target transcription regulatory network with the 42 key genes in HCC by Cytoscape 3.6.0 software (Figure4). The TF-target network complex contained 42 key genes and 7 TFs. We performed enrichment analysis for the 42 key genes; the results were followed in Figure 6.Figure 4
TF-related regulatory network. (a) TF-target network of these 42 key genes in HCC. (b) The brief framework of this network. TF-target network consisted of 42 nodes and 62 edges. The ellipses in the TF-gene network represented mRNAs with red (upregulated) and green (downregulated), and the diamonds represented TFs. The ellipses with a number were the clustered genes in the brief framework and the number of genes is shown inside. The interaction relationship between TFs and mRNAs were represented by arrows, and the direction of the arrow was from the source to the target. Different colors in the lines represented the different interaction relationship between the TFs and targets: red was for activation, green for repression, and grey for unknown.
(a) (b)Subsequently, the HCC and adjacent tissue gene expression profile of GSE64632 were analyzed by GEO2R. We use p<0.05 and |logFC|>1 as cut-off criterion to obtain 703 DEMs, including 452 upregulated DEMs and 251 downregulated DEMs. The interactions between TFs and microRNAs were predicted by miRTarBase. And we performed intersection analysis between the miRNA-targeting genes and DEMs in the microRNA assay data from GEO. Based on the above TRRUST version 2 database and miRTarBase, we constructed miRNA-TF-mRNA regulatory network by Cytoscape software (Figure5). The results showed that some genes were coregulated by the same miRNA or a few miRNAs, which suggested that miRNAs could play an important role in the progression of HCC. Using these coregulatory genes combined with regulatory networks, we constructed the miRNA-TF-mRNA regulatory network, which was shown in Figure 5. In addition, The Prediction results of differentially expressed miRNA targets revealed that EGR1 and FOS are coregulated by 2 microRNAs including hsa-miR-181a-5p and hsa-miR-192-5p, which were upregulated in HCC based on microRNA assay data. In the targeting relationship between miR-181a-5p and EGR1, FOS has been experimentally validated, especially.Figure 5
miRNA-TF-mRNA regulatory network for HCC. The squares in the network represented miRNAs, and the ellipses represented mRNAs, and the diamonds represented TFs. The nodes in red were upregulated, whereas the nodes in green were downregulated. The interaction relationship among TFs, mRNAs, and miRNAs was represented by arrows, and the direction of the arrow was from the source to the target. Different colors in the lines represented the different interaction relationship among TFs, mRNAs, and miRNAs: red was for activation, green for repression, and grey for unknown.Figure 6
Enrichment analysis of key genes.


### 3.3. Functional and Pathway Analysis of the Key Genes

In the current study, we performed enrichment analysis for the 42 key genes which were differentially expressed in HCC tissues. And the results were shown in Figure6 and Table 7. And the most enriched KEGG pathway of key genes was p53 pathway.Table 7
Enrichment analysis of key genes.
(a)
Top 10 enriched biological processes of key genes
GO_Name GO_ID P value Count response to chemical GO:0042221 4.35E-16 29 cellular response to chemical stimulus GO:0070887 8.93E-16 25 response to oxygen-containing compound GO:1901700 4.92E-14 18 response to drug GO:0042493 7.49E-14 12 response to stimulus GO:0050896 2.07E-13 38 response to organic substance GO:0010033 2.53E-13 23 negative regulation of cellular process GO:0048523 6.20E-13 26 regeneration GO:0031099 1.64E-12 9 animal organ development GO:0048513 1.78E-12 22 negative regulation of biological process GO:0048519 6.23E-12 26  (b)
Top 10 enriched molecular functions of key genes
GO_Name GO_ID P value Count protein binding GO:0005515 4.10E-08 32 oxidoreductase activity, acting on CH or CH2 groups, quinone or similar compound as acceptor GO:0033695 5.81E-06 2 caffeine oxidase activity GO:0034875 5.81E-06 2 cyclin-dependent protein kinase activity GO:0097472 8.20E-06 3 growth factor binding GO:0019838 3.10E-05 4 extracellular matrix binding GO:0050840 4.22E-05 3 chemokine activity GO:0008009 4.81E-05 3 oxidoreductase activity, acting on CH or CH2 groups GO:0016725 5.31E-05 2 receptor binding GO:0005102 5.38E-05 10 oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen GO:0016712 6.90E-05 3  (c)
Top 10 enriched cell components of key genes
GO_Name GO_ID P value Count extracellular region part GO:0044421 1.59E-09 20 extracellular region GO:0005576 1.97E-09 22 extracellular space GO:0005615 2.94E-09 14 membrane-bounded organelle GO:0043227 1.77E-08 38 cytoplasm GO:0005737 2.42E-08 36 platelet alpha granule lumen GO:0031093 5.55E-07 4 organelle GO:0043226 1.09E-06 38 platelet alpha granule GO:0031091 1.90E-06 4 cytoplasmic part GO:0044444 2.72E-06 30 membrane-bounded vesicle GO:0031988 3.14E-06 15  (d)
Top 10 enriched KEGG pathway of key genes
Pathway Name Pathway ID P value Count p53 signaling pathway hsa04115 8.24E-09 7 Cell cycle hsa04110 1.31E-04 5 Retinol metabolism hsa00830 1.31E-04 4 AMPK signaling pathway hsa04152 1.36E-04 5 Drug metabolism - cytochrome P450 hsa00982 1.65E-04 4 Metabolism of xenobiotics by cytochrome P450 hsa00980 2.06E-04 4 Progesterone-mediated oocyte maturation hsa04914 6.36E-04 4 Oocyte meiosis hsa04114 1.49E-03 4 Steroid hormone biosynthesis hsa00140 1.65E-03 3 FoxO signaling pathway hsa04068 2.04E-03 4Table 8
Primer sequences used for qRT-PCR amplification.
Primer 5′ > 3′ EGR1 CACCTGACCGCAGAGTCTTT GAGTGGTTTGGCTGGGGTAA  FOS GGGGCAAGGTGGAACAGTTA AGTTGGTCTGTCTCCGCTTG  FOSB GCGCCGGGAACGAAATAAAC AACTGATCTGTCTCCGCCTG  CCNB1 AATAAGGCGAAGATCAACATGGC TTTGTTACCAATGTCCCCAAGAG  CCNB2 CCGACGGTGTCCAGTGATTT TGTTGTTTTGGTGGGTTGAACT  CHEK1 ATATGAAGCGTGCCGTAGACT TGCCTATGTCTGGCTCTATTCTG  GAPDH GACAGTCAGCCGCATCTTCT ACCAAATCCGTTGACTCCGA  hsa-U6 CTCGCTTCGGCAGCACA AACGCTTCACGAATTTGCGT  hsa-miR-181a-5p ACGCTGACCTATGAATTGACAGCCWe queried all the samples from TCGA liver hepatocellular carcinoma (LIHC) with RNA-seq v2 data in our study. The boxplots showed that the expression level of EGR1, FOS, and FOSB was significantly lower in primary tumor than that in the normal liver for LIHC patients from TCGA (p<0.001) (Figure7(a)). The overall survival analysis of the three TFs was performed using the Kaplan-Meier curve. We found that the overexpression of FOSB in HCC was associated with reduced survival in cases with three TFs high expression, compared with the remaining cases with low/medium expression (Figure 7(b)). Through correlation analysis, we found that there was a strong positive correlation between the expression level of EGR1 and FOS mRNA. And there existed strong positive correlation toward coexpression of FOS and FOSB among target samples. The correlation between the expression level of EGR1 and FOSB mRNA was moderate (Figure 7(c)).Figure 7
Study of the clinical association of EGR1, FOS, FOSB, CCNB1, CCNB2, and CHEK1 with the clinicopathologic parameters of hepatocellular carcinoma. (a) Boxplots depicting RNA expression levels of key genes in HCC (n = 371) versus nonmalignant liver (n = 50) from TCGA. (b) Kaplan-Meier plots comparing the overall survival rates in HCC cases (n=365) with high expression or without low/medium expression. The data was recruited from UALCAN. P<0.05 was considered statistically significant. (c) Correlation analysis of three TFs and three p53 markers. The data was recruited from Linkedomics (http://www.linkedomics.org).
(a) (b) (c)Moreover, CCNB1, CCNB2, and CHEK1, which were three of key genes in the TF-target network, were involved in the p53 pathway. The boxplots showed that the expression levels of CCNB1, CCNB2, and CHEK1 were significantly higher in primary tumor than those in the normal liver for LIHC patients from TCGA (p<0.001). The overall survival rates of patients with high expression of CCNB1, CCNB2, and CHEK1 were all significantly lower than those of patients with low/medium expression. The correlation analysis results showed that there existed strong positive correlation among CCNB1, CCNB2, and CHEK1 mRNA expression (Figure7(c)).


### 3.4. qRT-PCR

We performed qRT-PCR to examine the expression of three differentially expressed TFs including EGR1, FOS, and FOSB mRNA in HCC. The relative expressions of EGR1, FOS, and FOSB mRNA were 0.493±0.558-, 0.494±0.476-, and 0.500±0.551-fold downregulated in 20 tumor tissues versus adjacent nontumor tissues, respectively (Figure8). Similarly, we performed qRT-PCR to explore the relationship among miR-181a-5p, TFs, and the key markers of p53 signaling pathways including CCNB1, CCNB2, and CHEK1. The relative expressions of CCNB1, CCNB2, and CHEK1 mRNA were 3.938±3.887-, 3.225±3.388-, and 3.186±3.508-fold upregulated in 20 tumor tissues versus adjacent nontumor tissues, respectively (Figure 8). And the relative expression of hsa-miR-181a-5p was 1.694±1.236-fold upregulated in 20 tumor tissues versus adjacent nontumor tissues.Figure 8
Validation of key genes and hsa-miR-181a-5p expression in 20 pairs of HCC and adjacent nontumor tissues by qRT-PCR. Detection of EGR1, FOS, FOSB, CCNB1, CCNB2, and CHEK1 mRNA expression and hsa-miR-181a-5p expression in HCC versus adjacent nontumor tissues was performed using qRT-PCR. Levels of EGR1, FOS, and FOSB mRNA were 0.493±0.558-, 0.494±0.476-, and 0.500±0.551-fold downregulated in tumor tissues, respectively, compared to those in the adjacent nontumor ones. And the levels of CCNB1, CCNB2, and CHEK1 mRNA were 3.938±3.887-, 3.225±3.388-, and 3.186±3.508-fold upregulated. The relative expression of hsa-miR-181a-5p was 1.694±1.236-fold upregulated.∗p<0.05, ∗∗p<0.01, and ∗∗∗p<0.001.


## 3.1. Identification and Enrichment Analysis of DEGs in HCC

HCC and adjacent nontumor tissue gene expression profile of GSE84402, GSE33006, GSE84005, and GSE12941 were obtained from NCBI-GEO. We used p<0.05 and |logFC|>1 as cut-off criterion. After integrated bioinformatical analysis, a total of 237 differentially expressed genes were identified from the four profile datasets, including 57 upregulated genes and 180 downregulated genes (Tables1–3, Figure 2).Table 1
The information of expression profiles.
Series Platform  Tissues Adjacent tissues Tumor tissues Reference GSE84402 GPL570 mRNA Liver 14 14 [4] GSE33006 GPL570 mRNA Liver 3 3 [18] GSE84005 GPL5175 mRNA Liver 38 38 - GSE12941 GPL5175 mRNA Liver 10 10 [19] GSE64632 GPL18116 miRNA Liver 3 3 [20]Table 2
The differentially expressed genes of downloaded expression profiles (P value<0.05, FC>2).
Upregulation Downregulation Total GSE84402 2232 1363 3595 GSE33006 1680 2696 4376 GSE84005 481 530 1011 GSE12941 175 385 560 Co-overexpression DEGs 57 180 237Table 3
The common DEGs of four gene expression profiles (P value<0.05, FC>2).
Common DEGs Gene symbol Downregulated DEGs CLEC4M, CXCL14, NAT2, FOSB, CFP, CRHBP, MARCO, CYP1A2, CLEC4G, FCN3, CNDP1, ABCA8, APOF, GYS2, GLS2, BMPER, ADRA1A, OIT3, VIPR1, GBA3, CDA, SDS, MT1G, SHBG, CYP3A4, HAMP, HAO2, GREM2, CYP39A1, HGFAC, C7, ANGPTL1, TIMD4, DPT, FOS, INMT, NNMT, TAT, C9, DNASE1L3, AKR1D1, IGJ, CYP8B1, PGLYRP2, TDO2, GPM6A, HHIP, BBOX1, CYP2B6, MT1X, RDH16, NPY1R, ZG16, DCN, GNMT, F9, C8A, HPGD, CETP, MFAP4, PLAC8, SRD5A2, AFM, ITIH4, LIFR, STAB2, HGF, C14ORF68, SRPX, TMEM27, ASPA, KLKB1, EGR1, ALDH8A1, SAA4, PHYHD1, CYP2C9, C6, DIRAS3, GHR, GNA14, ALPL, CD5L, PCK1, STEAP4, BMP5, IGF1, AGXT2, PZP, PON3, ACMSD, EPHA3, HPX, PRG4, RBMS3, ST3GAL6, FBP1, HBB, HAL, FTCD, ENO3, COLEC10, STEAP3, PROZ, CXCL12, SPP2, PDE7B, P2RY13, SLC38A4, ECM1, C8ORF4, SLC7A2, PRELP, C8B, CCL19, PTPRB, ITGA9, CHST4,ENG, RAMP3, GDF2, GPR128, ABCA9, CCDC3, VSIG4, GZMK, RND3, MBL2, LAMA2, EMR1, SOCS2, F8, CCBE1, CD69, CD163, PTGIS, FXYD1, DNAJC12, CA2, UROC1, GCH1, WDR72, RASGEF1B, COLEC11, IGFBP3, RGN, IFITM1, FBLN5, CYP2C18, CD4, GCDH, PKHD1, OLFML3, FGFR2, TEK, CPS1, TRIM22, ID1, GNAO1, SLC5A1, DHODH, PLSCR4, THBS1, ADH6, VNN1, HSD17B2, AXL, ZFP36, BGN, PDCD1LG2, SLC17A2, ASPN, SERPINA4, EHD3, FXYD6, IGFBP1, ACADSB, ADAMTS1, SFRP5, C1ORF162 Upregulated DEGs ENAH, TXNRD1, MCM4, GNPAT, FAM38B, SORT1, SRXN1, ECT2, TARBP1, RGS5, RFX5, TKT, RAP2A, TP53I3, TUBG1, ITGA6, PRKDC, AKR1C3, CCT3, TLCD1, ATAD2, HSPB1, CENPF, NUSAP1, SLC44A5, SPP1, CHEK1, CDH13, PEG10, MKI67, ASPM, ROBO1, CCNA2, BUB1, CKAP2, CDC6, UCK2, DTL, AURKA, TPX2, MDK, DEPDC1, GPC3, CAP2, UBE2T, PTTG1, CDKN3, RRM2, PRC1, FOXM1, CCNB2, SPINK1, ANLN, KIF4A, PBK, CCNB1, GINS1Figure 2
Identification of 237 common DEGs from the four cohort profile data sets. (a) 57 commonly upregulated DEGs. (b) 180 commonly downregulated DEGs. Different color areas represented different datasets. The cross areas meant the commonly changed DEGs. DEGs were identified with classical t-test; statistically significant DEGs were defined with p<0.05 and |logFC|>1 as the cut-off criterion.
(a) (b)The enrichment analysis of DEGs was performed. As shown in Table4 and Figure 3, in the biological process group, the DEGs were mainly enriched in response to chemical. In the cellular component group, the DEGs were most enriched in the extracellular region. In the molecular function group, the DEGs were most enriched in protein binding. And the most enriched KEGG pathway was complement and coagulation cascades.Table 4
Enrichment analysis of DEGs.
(a)
Top 10 enriched biological processes of DEGs
GO_Name GO_ID P value Count response to chemical GO:0042221 8.60E-31 103 cellular response to chemical stimulus GO:0070887 9.67E-27 80 response to stimulus GO:0050896 2.83E-25 158 single-organism metabolic process GO:0044710 1.07E-24 129 response to organic substance GO:0010033 1.49E-22 75 small molecule metabolic process GO:0044281 2.85E-22 76 single-multicellular organism process GO:0044707 6.06E-20 109 negative regulation of biological process GO:0048519 6.44E-20 91 negative regulation of cellular process GO:0048523 8.35E-20 86 regulation of response to stimulus GO:0048583 9.64E-20 82  (b)
Top 10 enriched cell components of DEGs
GO_Name GO_ID P value Count extracellular region GO:0005576 1.17E-37 110 extracellular region part GO:0044421 3.22E-35 96 membrane-bounded vesicle GO:0031988 3.92E-29 83 extracellular space GO:0005615 5.47E-29 60 vesicle GO:0031982 5.72E-29 89 extracellular vesicle GO:1903561 7.65E-27 67 extracellular organelle GO:0043230 7.96E-27 67 extracellular membrane-bounded organelle GO:0065010 2.19E-26 66 extracellular exosome GO:0070062 3.55E-26 66 cytoplasmic part GO:0044444 6.80E-14 136  (c)
Top 10 enriched molecular functions of DEGs
GO_Name GO_ID P value Count protein binding GO:0005515 4.15E-21 146 identical protein binding GO:0042802 2.61E-09 30 protein homodimerization activity GO:0042803 3.74E-09 22 oxygen binding GO:0019825 1.97E-08 8 growth factor binding GO:0019838 4.52E-08 10 extracellular matrix binding GO:0050840 6.88E-08 7 oxidoreductase activity GO:0016491 7.54E-08 32 binding GO:0005488 1.17E-07 186 protein dimerization activity GO:0046983 1.97E-07 27 oxidoreductase activity, acting on CH or CH2 groups GO:0016725 2.42E-07 4  (d)
Top 10 enriched KEGG pathways of DEGs
Pathway Name Pathway ID P value Count Complement and coagulation cascades hsa04610 1.79E-06 10 Metabolic pathways hsa01100 1.78E-05 43 p53 signaling pathway hsa04115 3.88E-05 8 Prion diseases hsa05020 3.97E-05 6 Histidine metabolism hsa00340 6.78E-05 5 Retinol metabolism hsa00830 1.93E-04 7 Steroid hormone biosynthesis hsa00140 6.97E-04 6 Tryptophan metabolism hsa00380 8.26E-04 5 Cell cycle hsa04110 2.17E-03 8 Caffeine metabolism hsa00232 3.35E-03 2Figure 3
Enrichment analysis of DEGs.


## 3.2. Identification and Functional Analysis of Key Differentially Expressed Genes, Construction of TF-Related Networks in HCC

We intersected 237 common DEGs and curated 36 cancer-related TF families to get three differentially expressed TFs (EGR1, FOS, and FOSB) in HCC (Table5). We searched the target genes and regulators of the three TFs through TRRUST version 2 database. Based on the results, we found that STAT6, a cancer-related TF, regulated FOS and FOSB, while FOS and EGR1 were coregulated by 2 cancer-related TFs, including BRCA1 and SP1, respectively. Similarly, TP53 and NKFB1 were coregulated by FOS and EGR1, respectively. Since FOS and EGR1 are linked to each other by 4 cancer-related TFs including BRCA1, SP1, TP53, and NKFB1, we also included these four TFs in the subsequent data of network construction. We screened the 4 TFs in the TRRUST version 2 database query to get each TF target genes and the DEGs from four profile datasets in GEO and TF-target genes for intersection analysis to obtain 42 differentially expressed target genes of the TFs (key genes) (Table 6), which laid the foundation for our next step to construct the gene signaling regulatory network in HCC.Table 5
Curated 36 cancer-related TF families.
Family Full Name Members (Official Gene Symbols) AP1 Activator Protein 1 FOS, FOSB, JUN, JUNB, JUND  AP2 Activator Protein 2 TFAP2A, TFAP2B, TFAP2C, TFAP2D, TFAP2E  AR Androgen Receptor AR  ATF Activating Transcription Factor ATF1 – 7  BCL B-cell CLL/lymphoma BCL3, BCL6  BRCA breast cancer susceptibility protein BRCA1 – 3  CEBP CCAAT/enhancer binding protein CEBPA, CEBPB, CEBPD, CEBPE, CEBPG  CREB cAMP responsive element binding protein CREB1 - 5, CREM  E2F E2F transcription factor E2F1 – 7  EGR early growth response protein EGR1 – 4  ELK member of ETS oncogene family ELK1, ELK3, ELK4  ER Estrogen Receptor ESR1, ESR2  ERG ets-related gene ERG  ETS ETS-domain transcription factor ETS1, ETS2, ETV4, SPI1  FLI1 friend leukemia integration site1 FLI1  GLI glioma-associated oncogene homolog GLI1 - 4  HIF Hypoxia-inducible factor HIF1A, ARNT, EPAS1, HIF3A  HLF hepatic leukemia factor HLF  HOX homeobox gene HOXA, HOXB, HOXD series, CHX10, MSX1, MSX2, TLX1, PBX2  LEF lymphoid enhancing factor LEF1  MYB myeloblastosis oncogene MYB, MYBL1, MYBL2  MYC myelocytomatosis viral oncogene homolog MYC  NFI nuclear factor I; CCAAT-binding transcription factor NFIA, NFIB, NFIC, NFIX  NFKB Nuclear factor kappa B, reticuloendotheliosis oncogene NFKB1, NFKB2, RELA, RELB, REL  OCT Octamer binding proteins POU2F1 - 3, POU3F1 - 2, POU5F1  p53 P53 family TP53, TP73L, TP73  PAX paired box gene PAX1 – 9  PPAR Peroxisome proliferator-activated receptor PPARA, PPARD, PPARG  PR Progesterone Receptor PGR  RAR retinoic acid receptor RARA, RARB, RARG  SMAD Mothers Against Decapentaplegic homolog SMAD1 - 9  SP sequence-specific transcription factor SP1 – 8  STAT signal transducer and activator of transcription STAT1 - 6  TAL1 T-cell acute lymphocytic leukemia-1 protein TAL1  USF upstream stimulatory factor USF1, USF2  WT1 Wilms tumor 1 (zinc finger protein) WT1Table 6
Differentially expressed target genes of the transcription factor.
Gene symbol differentially expressed target genes of the transcription factor CXCL14, FOSB, CYP1A2, ADRA1A, CYP3A4, FOS, PGLYRP2, CYP2B6, CETP, HGF, EGR1, PCK1, IGF1, FBP1, CXCL12, CCL19, CHST4, ENG, F8, CD69, CD163, IGFBP3, RGN, TRIM22, ID1, THBS1, ADH6, HSD17B2, BGN, ASPN, IGFBP1, CCNB1, TP53I3, PRC1, MKI67, ASPM, FOXM1, PTTG1, CHEK1, CCNB2, CCNA2, SPP1Based on the above TRRUST version 2 database, we constructed the TF-target transcription regulatory network with the 42 key genes in HCC by Cytoscape 3.6.0 software (Figure4). The TF-target network complex contained 42 key genes and 7 TFs. We performed enrichment analysis for the 42 key genes; the results were followed in Figure 6.Figure 4
TF-related regulatory network. (a) TF-target network of these 42 key genes in HCC. (b) The brief framework of this network. TF-target network consisted of 42 nodes and 62 edges. The ellipses in the TF-gene network represented mRNAs with red (upregulated) and green (downregulated), and the diamonds represented TFs. The ellipses with a number were the clustered genes in the brief framework and the number of genes is shown inside. The interaction relationship between TFs and mRNAs were represented by arrows, and the direction of the arrow was from the source to the target. Different colors in the lines represented the different interaction relationship between the TFs and targets: red was for activation, green for repression, and grey for unknown.
(a) (b)Subsequently, the HCC and adjacent tissue gene expression profile of GSE64632 were analyzed by GEO2R. We use p<0.05 and |logFC|>1 as cut-off criterion to obtain 703 DEMs, including 452 upregulated DEMs and 251 downregulated DEMs. The interactions between TFs and microRNAs were predicted by miRTarBase. And we performed intersection analysis between the miRNA-targeting genes and DEMs in the microRNA assay data from GEO. Based on the above TRRUST version 2 database and miRTarBase, we constructed miRNA-TF-mRNA regulatory network by Cytoscape software (Figure5). The results showed that some genes were coregulated by the same miRNA or a few miRNAs, which suggested that miRNAs could play an important role in the progression of HCC. Using these coregulatory genes combined with regulatory networks, we constructed the miRNA-TF-mRNA regulatory network, which was shown in Figure 5. In addition, The Prediction results of differentially expressed miRNA targets revealed that EGR1 and FOS are coregulated by 2 microRNAs including hsa-miR-181a-5p and hsa-miR-192-5p, which were upregulated in HCC based on microRNA assay data. In the targeting relationship between miR-181a-5p and EGR1, FOS has been experimentally validated, especially.Figure 5
miRNA-TF-mRNA regulatory network for HCC. The squares in the network represented miRNAs, and the ellipses represented mRNAs, and the diamonds represented TFs. The nodes in red were upregulated, whereas the nodes in green were downregulated. The interaction relationship among TFs, mRNAs, and miRNAs was represented by arrows, and the direction of the arrow was from the source to the target. Different colors in the lines represented the different interaction relationship among TFs, mRNAs, and miRNAs: red was for activation, green for repression, and grey for unknown.Figure 6
Enrichment analysis of key genes.


## 3.3. Functional and Pathway Analysis of the Key Genes

In the current study, we performed enrichment analysis for the 42 key genes which were differentially expressed in HCC tissues. And the results were shown in Figure6 and Table 7. And the most enriched KEGG pathway of key genes was p53 pathway.Table 7
Enrichment analysis of key genes.
(a)
Top 10 enriched biological processes of key genes
GO_Name GO_ID P value Count response to chemical GO:0042221 4.35E-16 29 cellular response to chemical stimulus GO:0070887 8.93E-16 25 response to oxygen-containing compound GO:1901700 4.92E-14 18 response to drug GO:0042493 7.49E-14 12 response to stimulus GO:0050896 2.07E-13 38 response to organic substance GO:0010033 2.53E-13 23 negative regulation of cellular process GO:0048523 6.20E-13 26 regeneration GO:0031099 1.64E-12 9 animal organ development GO:0048513 1.78E-12 22 negative regulation of biological process GO:0048519 6.23E-12 26  (b)
Top 10 enriched molecular functions of key genes
GO_Name GO_ID P value Count protein binding GO:0005515 4.10E-08 32 oxidoreductase activity, acting on CH or CH2 groups, quinone or similar compound as acceptor GO:0033695 5.81E-06 2 caffeine oxidase activity GO:0034875 5.81E-06 2 cyclin-dependent protein kinase activity GO:0097472 8.20E-06 3 growth factor binding GO:0019838 3.10E-05 4 extracellular matrix binding GO:0050840 4.22E-05 3 chemokine activity GO:0008009 4.81E-05 3 oxidoreductase activity, acting on CH or CH2 groups GO:0016725 5.31E-05 2 receptor binding GO:0005102 5.38E-05 10 oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen GO:0016712 6.90E-05 3  (c)
Top 10 enriched cell components of key genes
GO_Name GO_ID P value Count extracellular region part GO:0044421 1.59E-09 20 extracellular region GO:0005576 1.97E-09 22 extracellular space GO:0005615 2.94E-09 14 membrane-bounded organelle GO:0043227 1.77E-08 38 cytoplasm GO:0005737 2.42E-08 36 platelet alpha granule lumen GO:0031093 5.55E-07 4 organelle GO:0043226 1.09E-06 38 platelet alpha granule GO:0031091 1.90E-06 4 cytoplasmic part GO:0044444 2.72E-06 30 membrane-bounded vesicle GO:0031988 3.14E-06 15  (d)
Top 10 enriched KEGG pathway of key genes
Pathway Name Pathway ID P value Count p53 signaling pathway hsa04115 8.24E-09 7 Cell cycle hsa04110 1.31E-04 5 Retinol metabolism hsa00830 1.31E-04 4 AMPK signaling pathway hsa04152 1.36E-04 5 Drug metabolism - cytochrome P450 hsa00982 1.65E-04 4 Metabolism of xenobiotics by cytochrome P450 hsa00980 2.06E-04 4 Progesterone-mediated oocyte maturation hsa04914 6.36E-04 4 Oocyte meiosis hsa04114 1.49E-03 4 Steroid hormone biosynthesis hsa00140 1.65E-03 3 FoxO signaling pathway hsa04068 2.04E-03 4Table 8
Primer sequences used for qRT-PCR amplification.
Primer 5′ > 3′ EGR1 CACCTGACCGCAGAGTCTTT GAGTGGTTTGGCTGGGGTAA  FOS GGGGCAAGGTGGAACAGTTA AGTTGGTCTGTCTCCGCTTG  FOSB GCGCCGGGAACGAAATAAAC AACTGATCTGTCTCCGCCTG  CCNB1 AATAAGGCGAAGATCAACATGGC TTTGTTACCAATGTCCCCAAGAG  CCNB2 CCGACGGTGTCCAGTGATTT TGTTGTTTTGGTGGGTTGAACT  CHEK1 ATATGAAGCGTGCCGTAGACT TGCCTATGTCTGGCTCTATTCTG  GAPDH GACAGTCAGCCGCATCTTCT ACCAAATCCGTTGACTCCGA  hsa-U6 CTCGCTTCGGCAGCACA AACGCTTCACGAATTTGCGT  hsa-miR-181a-5p ACGCTGACCTATGAATTGACAGCCWe queried all the samples from TCGA liver hepatocellular carcinoma (LIHC) with RNA-seq v2 data in our study. The boxplots showed that the expression level of EGR1, FOS, and FOSB was significantly lower in primary tumor than that in the normal liver for LIHC patients from TCGA (p<0.001) (Figure7(a)). The overall survival analysis of the three TFs was performed using the Kaplan-Meier curve. We found that the overexpression of FOSB in HCC was associated with reduced survival in cases with three TFs high expression, compared with the remaining cases with low/medium expression (Figure 7(b)). Through correlation analysis, we found that there was a strong positive correlation between the expression level of EGR1 and FOS mRNA. And there existed strong positive correlation toward coexpression of FOS and FOSB among target samples. The correlation between the expression level of EGR1 and FOSB mRNA was moderate (Figure 7(c)).Figure 7
Study of the clinical association of EGR1, FOS, FOSB, CCNB1, CCNB2, and CHEK1 with the clinicopathologic parameters of hepatocellular carcinoma. (a) Boxplots depicting RNA expression levels of key genes in HCC (n = 371) versus nonmalignant liver (n = 50) from TCGA. (b) Kaplan-Meier plots comparing the overall survival rates in HCC cases (n=365) with high expression or without low/medium expression. The data was recruited from UALCAN. P<0.05 was considered statistically significant. (c) Correlation analysis of three TFs and three p53 markers. The data was recruited from Linkedomics (http://www.linkedomics.org).
(a) (b) (c)Moreover, CCNB1, CCNB2, and CHEK1, which were three of key genes in the TF-target network, were involved in the p53 pathway. The boxplots showed that the expression levels of CCNB1, CCNB2, and CHEK1 were significantly higher in primary tumor than those in the normal liver for LIHC patients from TCGA (p<0.001). The overall survival rates of patients with high expression of CCNB1, CCNB2, and CHEK1 were all significantly lower than those of patients with low/medium expression. The correlation analysis results showed that there existed strong positive correlation among CCNB1, CCNB2, and CHEK1 mRNA expression (Figure7(c)).


## 3.4. qRT-PCR

We performed qRT-PCR to examine the expression of three differentially expressed TFs including EGR1, FOS, and FOSB mRNA in HCC. The relative expressions of EGR1, FOS, and FOSB mRNA were 0.493±0.558-, 0.494±0.476-, and 0.500±0.551-fold downregulated in 20 tumor tissues versus adjacent nontumor tissues, respectively (Figure8). Similarly, we performed qRT-PCR to explore the relationship among miR-181a-5p, TFs, and the key markers of p53 signaling pathways including CCNB1, CCNB2, and CHEK1. The relative expressions of CCNB1, CCNB2, and CHEK1 mRNA were 3.938±3.887-, 3.225±3.388-, and 3.186±3.508-fold upregulated in 20 tumor tissues versus adjacent nontumor tissues, respectively (Figure 8). And the relative expression of hsa-miR-181a-5p was 1.694±1.236-fold upregulated in 20 tumor tissues versus adjacent nontumor tissues.Figure 8
Validation of key genes and hsa-miR-181a-5p expression in 20 pairs of HCC and adjacent nontumor tissues by qRT-PCR. Detection of EGR1, FOS, FOSB, CCNB1, CCNB2, and CHEK1 mRNA expression and hsa-miR-181a-5p expression in HCC versus adjacent nontumor tissues was performed using qRT-PCR. Levels of EGR1, FOS, and FOSB mRNA were 0.493±0.558-, 0.494±0.476-, and 0.500±0.551-fold downregulated in tumor tissues, respectively, compared to those in the adjacent nontumor ones. And the levels of CCNB1, CCNB2, and CHEK1 mRNA were 3.938±3.887-, 3.225±3.388-, and 3.186±3.508-fold upregulated. The relative expression of hsa-miR-181a-5p was 1.694±1.236-fold upregulated.∗p<0.05, ∗∗p<0.01, and ∗∗∗p<0.001.


## 4. Discussion

Although some progress has been made in the study of HCC, the exact mechanisms of occurrence and development in HCC are not yet clear. In the current study, we constructed networks related to transcription regulatory modes in HCC and performed functional and KEGG pathway analysis for the key genes. The enrichment analysis results showed that key genes in the TF regulatory network were enriched in the p53 signaling pathway. In addition, we performed qRT-PCR and verification in the TCGA Dataset to confirm the results based on bioinformatics analysis.FOS (also known as c-FOS) and FOSB were both protooncogenes, which belong to the activator protein 1(AP1) family [27]. Many studies have reported that FOS was involved in proliferation, migration, and invasion of some malignancies [28, 29]. In some cases, expression of the FOS gene has also been associated with apoptotic cell death [30]. However, there still existed some contradictions about the role that c-FOS might play in the tumor progression. Some studies showed that overexpression of FOS might be associated with inhibition of tumor [31, 32]. Oliveira-Ferrer L et al. reported that c-FOS overexpression increased the apoptotic potential of ovarian cancer cells and inhibited tumor growth and metastasis, which could be achieved by changing the adhesion of human ovarian cancer cell [33]. A study by Guo J et al. showed that the expression of c-FOS in the tumor tissues of pancreatic cancer seemed to be lower than that in adjacent nontumor tissues [34]. However, some other research showed that FOS could contribute to carcinogenesis [35, 36]. Therefore, the above research work indicated that c-FOS is expressed differently in different histological types and is closely related to the proliferation, differentiation, invasion, and apoptosis of tumor cells, which provided a new therapeutic target for cancer by regulating the expression of c-FOS. EGR1 (early growth response 1) belongs to EGR family of transcription factors that includes four members: EGR1, EGR2, EGR3, and EGR4. It is a nuclear protein and functions as a transcriptional regulator, which is a component of p53 signaling [37, 38]. Abundant studies found that expression of EGR1 was associated with HCC metastasis and proliferation [39, 40]. And a number of studies suggested that EGR1 exhibits prominent tumor-suppressive activity by activating major tumor suppressor factors, including transforming growth factor-β1, p53, p73, fibronectin, and PTEN [40, 41]. In the present study, the mRNA expression of FOS and EGR1 was lower in HCC tissues than that in nontumor adjacent tissues based on bioinformatics analysis and qPCR verification. And the low expression of FOS and EGR1 was negatively associated with overall survival of HCC patients based on TCGA data, although this correlation was not statistically significant. The results above suggested that EGR1, FOS, and FOSB, three cancer-related TFs, could play an important role in the progression of HCC.In the current study, we found that miR-181a-5p was upregulated in HCC based on miRNA microarray data and qRT-PCR verification. As a member of the miR-181 family, the level of miR-181a-5p was overexpressed in many cancers including breast cancers, multiple myeloma, pancreatic and gastric cancer, and hepatocellular carcinoma [42–45]. And some studies have reported that miR-181a was involved in the pathogenesis of HCC by inducing hepatocyte epithelial-mesenchymal transition and decreasing autophagy [46, 47]. A study by C. Zou et al. confirmed that miR-181a plays an important role in the progression of HCC autophagy-related modulator [48]. We found that FOS and EGR1 were both targeted by miR-181a-5p based on prediction results of key genes target, and the targeting relationship between miR-181a-5p and FOS, EGR1, has been confirmed in some studies [49, 50]. And based on correlation analysis from TCGA data, there was a tendency to a positive association between FOS and EGR1. Furthermore, the expression of FOS protein downregulated by miR-181a has been reported [51]. And another study revealed that aberrant EGR1 expression could be suppressed by miR-181a-5p directly [52].The role of the p53 pathway in the progression of HCC has been reported in the literature [53, 54]. As an important component of the p53 pathway, TP53 (tumor protein p53) is the common target of EGR1 and FOS. Studies found that increased EGR1 expression could activate p53 signaling pathway to induce apoptosis in HCC cells [55, 56]. Association between FOS and p53 signaling has also been reported, but the specific interaction relationship is not yet clear [57]. In the current study, we found that CCNB1, CCNB2, and CHEK1 were the target genes of TP53 based on bioinformatics analysis and literature confirmation. Cyclin B1 (CCNB1) and cyclin B2 (CCNB2) are important members of the cyclin family and are important cell cycle regulators related to G2/M detection sites [58, 59]. One of its important roles is to modulate and form a complex with cyclin-dependent kinase 1 (CDK1) to phosphorylate the substrate, initiate the cell to G2/M phase from G1/S, and promote mitosis [60]. Checkpoint kinase 1 (CHEK1), as a DNA damage sensor and cell death pathway stimulator, regulates the progression of the cell cycle at the S phase, G2/M checkpoint [61]. Damaged DNA activates CHEK1, causing cell cycle arrest and repairing damaged DNA; if extensive damage is not repaired, it induces apoptosis, thereby maintaining genome integrity and stability [62]. Due to its anti-injury effect, CHEK1 plays an important role in tumor development and apoptosis [63–65]. Many studies have reported that CCNB1, CCNB2, and CHEK1 were involved in p53 signaling pathway [66–68]. We found that mRNA levels of three p53 markers including CCNB1, CCNB2, and CHEK1 were significantly upregulated in 20 HCC tumor tissues versus nontumor adjacent tissues by bioinformatics analysis and experimental verification, and their mRNA expression levels were all negatively correlated with HCC patients’ survival rates (p<0.001). Therefore, we speculated that as an onco-microRNA, miR-181a-5p could inhibit expression of FOS and EGR1 to regulate p53 signaling pathway, which may be achieved through the upregulation of CCNB1, CCNB2, and CHEK1, thereby promoting the progression of HCC.


## 5. Conclusion

Using multiple cohorts profile datasets and integrated bioinformatical analysis, we identified commonly 237 DEGs, and finally found EGR1, FOS, and FOSB, 3 cancer-related TFs, which were downregulated in HCC. In addition, we constructed the transcription factor-related regulatory networks based on EGR1, FOS, and FOSB and identified possible miR-181a-5p→FOS/EGR1→p53 signaling pathway axis. It should be noted that this study examined the conclusion by qRT-PCR and bioinformatics analysis; further research needs to be done to explore more specific mechanisms. Notwithstanding its limitation, these findings significantly improved the understanding of underlying molecular mechanisms in HCC, and the key genes and pathways could be used as diagnostic and therapeutic targets and diagnostic biomarkers.



---

*Source: 1431396-2018-08-29.xml*