Identification of FOS as a Candidate Risk Gene for Liver Cancer by Integrated Bioinformatic Analysis

Liver cancer is a lethal disease that is associated with poor prognosis. In order to identify the functionally important genes associated with liver cancer that may reveal novel therapeutic avenues, we performed integrated analysis to profile miRNA and mRNA expression levels for liver tumors compared to normal samples in The Cancer Genome Atlas (TCGA) database. We identified 405 differentially expressed genes and 233 differentially expressed miRNAs in tumor samples compared with controls. In addition, we also performed the pathway analysis and found that mitogen-activated protein kinases (MAPKs) and G-protein coupled receptor (GPCR) pathway were two of the top significant pathway nodes dysregulated in liver cancer. Furthermore, by examining these signaling networks, we discovered that FOS (Fos proto-oncogene, AP-1 transcription factor subunit), LAMC2 (laminin subunit gamma 2), and CALML3 (calmodulin like 3) were the most significant gene nodes with high degrees involved in liver cancer. The expression and disease prediction accuracy of FOS, LAMC2, CALML3, and their interacting miRNAs were further performed using a HCC cohort. Finally, we investigated the prognostic significance of FOS in another HCC cohort. Patients with higher FOS expression displayed significantly shorter time to recurrence (TTR) and overall survival (OS) compared with patients with lower expression. Collectively, our study demonstrates that FOS is a potential prognostic marker for liver cancer that may reveal a novel therapeutic avenue in this lethal disease.


Introduction
Liver cancer is a lethal disease often caused by liver damage associated with virus infection, excessive alcohol, or other liver disease [1]. It ranks the sixth most common cancer worldwide and accounts for the fourth most common cause of cancer related death due to lack of effective therapies [2]. Currently, the most curative therapy for liver cancer is surgery, but most patients are not suitable for this treatment because of the advanced tumor stage at the time of diagnosis. Currently, there are limited targeted therapy agents or drug combinations that can extend the survival time of these patients [3].
In recent years, the studies of cancer genetics and genomic sequencing have contributed to deeper understanding of the underlying cause of liver cancer [4]. e gene expression profiling analysis has improved the classification of liver cancer subtypes [5]. e discovery of mutations in BRAF (B-Raf proto-oncogene, a serine/threonine kinase) and EGFR (epidermal growth factor receptor) by DNA sequencing contributed to the development of new targeted therapies for liver cancer patients carrying these mutations [6,7]. Based on gene expression patterns of liver cancer, MST1 (or so called "hepatocyte growth factor-like protein") was found to be lowly expressed in liver cancer sample induced by p53 [8], and its role in regulation of hepatocyte differentiation has been studied previously [9]. In addition, miR-122 was found to be expressed at lower levels in liver cancer cells compared to normal hepatocytes by gene expression profiling, thus contributing to its function on liver cancer cell migration and invasion [10]. miR-122 has been suggested as a diagnostic and prognostic marker for liver cancer. However, the molecule marker for liver cancer has not been fully investigated.
In our current study, we aimed to identify novel genes or pathways that may be important for liver cancer progression by performing integrated analyses of mRNA and miRNA in liver cancer patients. e global genomic analysis may be helpful to discover new cancer-associated genes and reveal the molecular basis important for liver cancer development.

Data Acquisition and Differential Expression Analysis.
e miRNA and mRNA expression datasets of liver cancer were downloaded from TCGA database (https://gdc-portal. nci.nih.gov/). ere were 372 liver cancer tissues and 50 tumor-adjacent tissues for miRNA expression profiles and 371 tumor samples and 50 adjacent tissues for mRNA expression dataset.
Compared with the tumor-adjacent tissues, the differentially expressed (DE) mRNAs (genes) and DE miRNAs in tumor tissues were identified by edgeR package in R. Corrected P value <0.05 and |log 2 FC(fold change)| >1.0 were considered as significant.

Risk Genes for Liver Cancer Analysis.
e targets of experiment-validated DE miRNAs were predicted based on miRecords, miRTarBase, and TarBase. e target genes were mapped to DE mRNAs, and the overlapped genes were defined as the risk genes of liver cancer.

Differentially Expressed Pathway
Analysis. KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www. genome.jp/kegg/pathway.htm) is a database associated with gene function annotation, which is comprised of four main databases including GENES database, PATHWAY database, LIGAND database, and BRITE database [11]. e PATH-WAY database is the collection of pathways for various biological processes. All the human pathways and their related genes were downloaded from the PATHWAY database. e pathway expression value in each sample was calculated based on the median expression level of each gene involved in pathway. e formula was listed as follows: where Path ik represents the expression level of pathway i in sample k and g 1−n represents the expression value of each gene in pathway i.
Compared with controls, the differentially expressed pathways with corrected P value <0.05 and |log 2 FC| ≥ 1 in tumor samples were analyzed by edge package in R. Bidirectional clustering analysis of the differentially expressed pathways was carried out by pheatmap package in R.

e Correlated Pathway Pairs.
e correlation between pathways in tumor and control samples was analyzed by Pearson correlation coefficient as follows: where P X,Y represents the expression correlation coefficient between pathway X and Y and X and Y represent the mean expression value of pathway X and Y, respectively. e pathway pairs with |correlation coefficient| >0.8 in both tumor and control samples were collected.

miRNA-Gene-Pathway Network Construction.
In order to obtain risk gene-pathway pairs, the common risk genes in both pathways of correlated pathway pair were collected. miRNA-gene-pathway network was constructed based on risk gene-pathway pair and miRNA-risk gene interactions. Subsequently, the properties of network topology, such as degree, average shortest path length, betweenness centrality, closeness centrality, clustering coefficient, and topological coefficient, were analyzed by "network analysis" function of Cytoscape software. e disease classification accuracy of significant gene and miRNAs was predicted by ten-fold cross validation.

2.7.
Patients and Follow-Up. 382 patients who were diagnosed with HCC and underwent surgery from January to December in 2011 in Zhongshan Hospital were recruited. e inclusion criteria were as follows: (a) no prior cancer treatment, (b) pathologically diagnosed HCC, with complete resection of all tumor nodules with margins confirmed free of cancer by histologic examination, and (c) availability of complete clinicopathologic and follow-up data. e Barcelona Clinic Liver Cancer (BCLC) staging system was used to assess tumor stage. Tumor differentiation was determined according to the Edmondson grading system. Approval for research protocol and use of human subjects was obtained from the research ethics committee of Zhongshan Hospital. Informed consent was obtained from each patient. Followup ended in August 2017. Time to recurrence (TTR) was defined as the interval between surgery and the diagnosis of any type of recurrence including intra or extrahepatic recurrence identified by MR or CT. OS was defined as the interval between treatment and death of any cause or the last observation date [12].

Tissue Microarray (TMA) and Immunohistochemistry (IHC).
Immunohistochemical staining was performed using the avidin-biotin-peroxidase complex method. Primary anti-human-FOS antibodies were added to the slides and incubated at 4°C overnight, followed by rehydration and microwave antigen retrieval. Subsequently, secondary antibody was incubated at 37°C for 30 minutes. Slides were stained with 3′3-diaminobenzidine tetra hydrochloride and counterstained with Mayer's hematoxylin. e assessment of immunohistochemical staining was performed by two independent pathologists, and discrepancies were resolved by consensus. e intensity of FOS staining was stratified as weak or strong.

Statistical Analysis.
Statistical analyses were performed using SPSS 20.0 software (IBM, Armonk, NY, USA). Experimental values for continuous variables are expressed as the mean ± standard error of the mean. Chi-squared tests, Fisher's exact probability tests, and Student's t-tests were used to evaluate the significance of differences between groups. If variances within groups were not homogeneous, a nonparametric Mann-Whitney test or Wilcoxon signedrank test was used. e relationships between FOS expression and TTR or OS were analyzed using Kaplan-Meier survival curves and log-rank tests, respectively. P < 0.05 was considered statistically significant [13].

Identification of Potentially Functionally Important mRNAs and miRNAs for Liver Cancer.
In order to examine the potentially important mRNAs and miRNAs in liver carcinogenesis, first we performed analysis of differentially expressed (DE) mRNAs or miRNAs from liver cancer patients compared to normal tissues in TCGA database by using corrected-P value < 0.05 and |log 2 FC| > 1.0; a total of 405 DE genes (273 upregulated ones and 132 downregulated ones) were identified based on gene expression profiling of liver tumors compared to normals. In addition, 233 DE miRNAs (39 upregulated miRNAs and 194 downregulated miRNAs) were obtained based on the miRNA expression profiling in the same dataset.
Next, to examine the functional network between miRNAs and their target genes, we obtained a total of 324475 experiment-validated miRNA-target interactions, among which there were 5559 DE miRNA-target interactions, including 86 DE miRNAs and 3378 target genes. After the target genes were mapped to DE genes, 35 DE miRNA-DE gene pairs were obtained, which included 22 DE miRNAs and 24 risk genes ( Figure 1). is is a reasonably focused miRNA and gene list that will be used for our subsequent analysis in liver cancer.

Examination of Potentially Important Pathways for Liver
Carcinogenesis. To determine the pathways dysregulated in liver cancer due to these differentially expressed mRNAs or miRNAs, we used the cutoff value of corrected P value < 0.05 and |log 2 FC| ≥ 1 to perform the gene set enrichment analysis with the KEGG pathway database. A total of 622 pathways were found to be significantly dysregulated in liver tumor samples compared to normal, which was presented as a heatmap plot in Supplementary Figure 1. erefore, these pathways might have significant roles in liver tumorigenesis.

Elucidation of miRNA-Gene-Pathway Network in Liver
Cancer. To examine the crosstalk between miRNA-mRNA dysregulated in liver cancer, which may yield more confidence on studying important oncogenic pathways in this cancer, we examined the miRNA-mRNA-pathway correlation. By using the |correlation coefficient| >0.8, we obtained 13424 correlated pathway pairs. After comparing the risk genes and pathway related genes, the miRNA-risk genepathway pair network was constructed, which was comprised of 1024 edges connecting 115 nodes ( Figure 2). en, the topological properties of nodes in network were achieved by Cytoscape software. MAPK signaling pathway 24   (degree � 53) and GPCR Pathway (degree � 51) were two of the most significant pathways in miRNA-gene-pathway network that showed most interactions with risk genes (Table 1). Further analyses of important genes in these pathways revealed that FOS (Fos proto-oncogene, AP-1 transcription factor subunit, degree � 54), LAMC2 (laminin subunit gamma 2, degree � 13), and CALML3 (calmodulin like 3, degree � 7) were the most significant gene nodes with high degrees indicating the pivotal role of these genes involved in liver cancer (Table 2). FOS was previously reported to be regulated by hsa-miR-221 and hsa-miR-222 and showed interactions with MAPK signaling pathway, JAK/STAT pathway, GPCR pathway, and insulin receptor pathway (Figure 3(a)). LAMC2 interacted with hsa-miR-199b and was  involved in the pathways such as MAPK signaling, ERK signaling, PI3K signaling, and PTEN pathway (Figure 3(b)). CALML3 interacted with hsa-miR-765 and played a regulatory role in GnRH signaling pathway, long-term potentiation, and Huntington's disease (Figure 3(c)).

Validation of Significant Genes and Interacting miRNAs in Liver Cancer.
To validate our miRNA-mRNA-pathway interaction network analysis, we examined the expression of FOS, LAMC2, CALML3, and their interacting miRNAs using 20 pairs of hepatocellular carcinoma samples and adjacent normal tissues by RT-PCR. Our results were highly consistent with our analysis from TCGA database (Figure 4(a)). Furthermore, the disease prediction accuracy of FOS, LAMC2, CALML3, and their interacting miRNAs was analyzed. As shown in Figure 4(b), FOS (0.9074), hsa-miR-199b (0.8772), and hsa-miR-221 (0.8222) displayed high receiver operating characteristic (ROC) values, among which FOS exhibited the highest prediction accuracy.

High FOS Expression Predicts Poor Prognosis in HCC
Patients. Based on our analysis above, we were motivated to study the relationship of FOS expression with the prognosis of HCC. To this end, a TMA containing 382 patients who underwent curative resection was immunostained with the previously characterized antibody against FOS. Basic pathological and clinical information for these 382 HCC patients enrolled is described in Table 3. ese HCC patients were divided into two groups according to the FOS expression level (high or low) determined by the IHC staining intensity (strong or weak) ( Figure 5(a)). Our statistical analysis showed that high expression of FOS was associated with Hepatitis B Virus (HBV) infection background, alphafetoprotein (AFP) level, and macrovascular invasion. Kaplan-Meier analysis revealed significantly shorter median overall survival (OS) in patients with higher FOS expression compared with those patients with lower FOS expression (24.7 months vs. not reached, * * P < 0.01, Figure 5(b)). Similarly, patients with higher FOS expression displayed significantly shorter time to recurrence (TTR) compared  Collectively, our results suggest that FOS can be used as an effective marker of prognosis in HCC.

Discussion
Tumor initiation is attributed to the process of genetic and epigenetic alterations of patients, which lead to gene expression alternations involved in tumor evolution. In this paper, we analyzed the gene expression and miRNA expression profiling of liver cancer by performing integrated bioinformatics analysis. e DE mRNAs and DE miRNAs were identified by comparing tumor tissues and controls followed by the miRNA-mRNA-pathway network construction.
e integrative analyses of miRNA-mRNA- pathway network for identifying functionally relevant genes in liver cancer progression are not frequently used. By performing these new analyses, we attempted to pinpoint driver genes that may contribute to the development of liver cancer. Our data showed that MAPK and GPCR pathways were two of the most significant nodes in miRNA-gene-pathway network. MAPKs are the family of serine-threonine kinases, which contain extracellular signal-regulated kinase (ERK), p38, and c-Jun NH2-terminal kinase (JNK). MAPK signaling pathway is activated by extracellular and intracellular stimuli and plays a key role in cell proliferation, differentiation, and other cellular activities [14]. MAPK signaling pathway has been reported to be implicated in pathogenesis of many different cancers. e dysregulation of MAPK signaling pathway in cancer development is frequently due to the mutations of signaling components involved in the pathway [15]. It is reported that the expression and activity of MAPK are significantly upregulated in primary liver cancer [16]. MAPK is found to be overexpressed in liver cancer patients and plays a key role in the growth and survival of liver cancer cells [17]. Targeting MAPK signaling pathway has been suggested to be one of the attractive candidates for cancer therapy. In the previous study, RAS/ MAPK signaling pathway activity is reported to be upregulated in medulloblastoma based on gene expression profiling [18]. MAPK signaling pathway has been found to be the most important signaling node based on lncRNA profiling in breast cancer [19]. Consistent with these findings in breast cancer, our study also suggests the significance of MAPK signaling in liver cancer.
GPCRs are the family of seven-transmembrane receptors and play key roles in signal transmission involved in various physiological functions. GPCRs are overexpressed in various cancers and play key roles in tumor cell proliferation. erefore, the dysfunction of GPCRs has been proven to BioMed Research International promote the progression and metastasis of cancers [20]. For example, GPCRs have been reported to be involved in prostate cancer development [21]. However, the direct evidence for the role of GPCR pathway in liver cancer is lacking. ere were some limited literatures suggesting that GPCRs are closely related to the activation of Hippo pathway [22] involved in liver growth and tumor formation [23]. All these literatures above suggest the significant roles of MAPK and GPCR pathways in liver cancer progression.
In our study, our data suggest that FOS is an important gene in signaling nodes involved in the interaction of MAPK and GPCR pathways. It is reported that c-FOS plays a key role in the signal transduction pathway which acts as a trans-activating as well as trans-repressing molecule [24]. e c-FOS is a mitogen responsive gene associated with cell proliferation. It is reported that c-FOS gene expression is activated by the combined effect of extracellular nucleotides and growth factors, which leads to increased calcium  levels and proliferation in breast cancer cells [24]. In addition, the stimulation of c-FOS gene expression is implicated in activation of ERK signaling. Furthermore, overexpression of c-FOS gene is mediated by the G proteincoupled receptor GPR30 through E2 (17β-estradiol) and phytoestrogens in breast cancer cells [25]. us, both ERK signaling and GPR30 signaling are implicated in cell proliferation of breast cancer induced by FOS gene. e significant role of FOS gene is also identified in other cancers, such as bladder cancer [26], lung cancer [27], and colon cancer [28]. However, the role of FOS gene in liver cancer remains largely unclear. In our present study, FOS gene exhibits the highest prediction accuracy of liver cancer among the significant nodes in network, which suggests the potential diagnostic and therapeutic implication of FOS in liver cancer.

Conclusion
In conclusion, the co-activation of MAPK and GPCR pathways may be involved in liver cancer progression. FOS may be the risk gene for liver cancer by playing important roles in MAPK and GPCR pathways. FOS could be considered to be a candidate gene for diagnosis and therapy for liver cancer.
Data Availability e data that support the findings of this study are available from the corresponding author upon reasonable request.