Deciphering the Novel Target Genes Involved in the Epigenetics of Hepatocellular Carcinoma Using Graph Theory Approach

Background Even after decades of research, cancer, by and large, remains a challenge and is one of the major causes of death worldwide. For a very long time, it was believed that cancer is simply an outcome of changes at the genetic level but today, it has become a well-established fact that both genetics and epigenetics work together resulting in the transformation of normal cells to cancerous cells. Objective In the present scenario, researchers are focusing on targeting epigenetic machinery. The main advantage of targeting epigenetic mechanisms is their reversibility. Thus, cells can be reprogrammed to their normal state. Graph theory is a powerful gift of mathematics which allows us to understand complex networks. Methodology In this study, graph theory was utilized for quantitative analysis of the epigenetic network of hepato-cellular carcinoma (HCC) and subsequently finding out the important vertices in the network thus obtained. Secondly, this network was utilized to locate novel targets for hepato-cellular carcinoma epigenetic therapy. Results The vertices represent the genes involved in the epigenetic mechanism of HCC. Topological parameters like clustering coefficient, eccentricity, degree, etc. have been evaluated for the assessment of the essentiality of the node in the epigenetic network. Conclusion The top ten novel epigenetic target genes involved in HCC reported in this study are cdk6, cdk4, cdkn2a, smad7, smad3, ccnd1, e2f1, sf3b1, ctnnb1, and tgfb1.


INTRODUCTION
C.H. Waddington's [1] definition of epigenetics illustrating it as the 'study of interactions between genes and the products thereof leading to the formation and existence of phenotype' was initially concerned with the role of epigenetics in fetal development. However, the definition of epigenetics has evolved over time as it is implicated in a wide variety of biological processes. The contemporary definition of epigenetics is 'the study of heritable and reversible changes in gene expression that occur independently of changes in the primary DNA sequence' i.e. somatic changes. Most of these heritable changes are established during differentiation and are stably maintained through multiple cycles of cell division, enabling cells to have distinct identities, while containing the same genetic information. This heritability of gene expression patterns is mediated by epigenetic modifications, which include methylation of CpG regions in DNA, post translational Histone modification, nucleosome positioning along with the DNA and micro-RNA modification.
Epigenome is an assembly of chemical entities that allocates and guides the genome to perform the intended function. DNA constitutes the human genome, whereas the *Address correspondence to this author at the Department of Bioinformatics and Applied Sciences, Indian Institute of Information Technology, Allahabad, Uttar Pradesh, India; Tel: 9236666060; E-mail: pritish@iiita.ac.in chemical compounds and proteins constitute the epigenome. These compounds and proteins attach to the DNA and perform functions like turning genes on or off as well as controlling the rate of production of proteins. Cancers are caused as a cumulative effect of the genome and epigenome [2,3].
For a very long time, it was assumed that genetics and epigenetics are two independent mechanisms [3,4]. Mutations in genetic and epigenetic mechanisms lead to cancer development and promote cancer progression [5]. DNA methylation, histone modification and micro-RNA changes are found to be the important biomarkers of the initial stages in many forms of cancer [6]. The fact that epigenetic aberrations, unlike genetic mutations, are potentially reversible and can be reprogrammed to their normal state by epigenetic therapy makes such initiatives promising and therapeutically relevant [7].
The three most studied epigenetic mechanisms include (a) DNA methylation, (b) histone modifications (both covalent and non-covalent) and (c) micro RNA or miRNA. The process of addition of methyl (CH 3 ) groups to the DNA molecule is known as DNA methylation. It can change DNA activity without altering the sequence. Hence, it is a crucial factor in epigenetics [8]. At normal levels, DNA methylation is required for normal development and inactivation of various processes like X-chromosome inactivation. When a CpG cluster region at the promoter site of a gene is methylated, expression of the gene is turned off.
Post-translational histone modifications are observed in many human cancers [9]. These modifications can be of many types. Some of the common Histone modifications include Acetylation, Methylation, Phosphorylation, Ubiquitylation and Sumoylation. miRNAs play very important cellular functions like regulation of mRNA activity [10]. They have proved to be early markers of cancers [11]. miR-NA control the rate of expression of various enzymes involved in epigenetics. Though the mechanism is not fully understood, the role of miRNA is an established fact.
In spite of its potential of early biomarker detection, epigenetic drug discovery poses a number of challenges ( Fig. 1) like less number of biological and chemical tools and assays for probing and screening of chemical compounds [12]. This needs to be addressed before epigenetic drugs make their way to the patients. . Challenges in drug discovery of epigenetic drugs. It includes the lack of biological reference compounds for assay development. Currently, there is a lack of knowledge about both the short term as well as the long term repercussions of the epigenetic therapies (toxicology). The need of the hour is new models and longer duration to study tumor biology in vivo. Presently, limited epigenetic proteins and antibodies of high quality (target selection) are available. (A higher resolution / colour version of this figure is available in the electronic copy of the article).
The mathematical discipline concerned with the study of complex networks is graph theory. A graph consists of points called vertices which are connected to each other via edges. This can be used to model many practical problems like biological systems [13]. Shifting from mainstream mathematics, it is now actively being used in fields like genomics, biological network analysis, electrical engineering, etc. The infinite scope of the combinatorial methods of graph theory leading to significant results in mathematics and other applications has been well-documented [14,15]. In the systems biological approach, the analysis of the interaction between various components is very important as it provides a deeper insight into the overall working of pathway. A centrality measure quantifies the importance of vertices by allotting ranks [16]. Various centrality measures are there like closeness centrality, degree centrality, motif-based centralities, etc. Depending upon the type of network, the type of centrality measure is chosen. Network analysis provides a powerful tool in order to understand the structure, function and the evolutionary patterns involved in the biological process [17].
Almost 90% of liver cancers result from Hepatocellular carcinoma (HCC) [18]. Approximately, 7.5 Lakhs of new cases of HCC per year occurs globally which makes it the 5 th common cause of cancers affecting human [18]. According to the report published by the International Agency of Research on Cancer (WHO), the male: female ratio for HCC in India is 4:1.
In this paper, we have merged the network of differentially expressed genes with the pathways most frequently affected in HCC to form an integrated network. This network was further validated, and topological parameters were calculated to quantify each node (Fig. 2). Fig. (2). The methodology implemented to identify novel targets.

Data Mining
The initial data for epigenetic network construction was incorporated from two resources.

GEO Datasets
Expression profiles of GSE18081 [19], GSE37988 [20], GSE44970 [21], GSE54503 [22] and GSE57956 [23] were obtained from the GEO database [24]. GPL570 platform was used for the analysis. GSE18081 data represents CpG site methylation of HCV-cirrhotic, HCV-HCC and normal liver tissues. GSE37988 represents the profile of methylation in Taiwanese HCC tumor and adjacent tissues. GSE44970 comprises genome-wide DNA methylation profiles of human hepatocellular carcinoma. GSE54503 consists of genomewide DNA methylation profiles altered in hepatocellular carcinoma using Infinium Human Methylation 450 Bead Chips. GSE57956 is the methylation profile of hepatocellular carcinoma. All the datasets were checked for median centricity to see if the datasets chosen are comparable or not.
These datasets were basically micro-array data of Homo sapiens. The cutoff value of adjusted p-value <0.05 and |log FC| value > 1.5 was chosen as the benchmark. The DEGs were then used for the construction of PPI network from STRING [25].

Literature Search
There are various epigenetic markers involved in cancer epigenetics. One of the most critical is DNA methylation. It could be both hypomethylation and hypermethylation. We searched for various pathways that are most frequently affected by HCC, incorporated genes and merged it with the DEG network. Some of the pathways are Wnt-beta catenin, p53, hedgehog, etc. The PPI network of the differentially expressed genes was obtained from the STRING database [25].

Formation of an Integrated Protein-protein Interaction (PPI) Network
Both the PPI network constructed from DEGs and the pathways searched from literature, were imported in Cytoscape 3.5.1 [26]. These networks were merged using the merge tool implemented in Cytoscape 3.5.1 using the union operation.

Network Validation
Epigenomic network was validated by comparing it to random network and examining if it followed the power-law degree distribution. A comparison with random network was done using the Network Randomizer 1.1.2 [27] plugin in Cytoscape 3.5.1. A modular approach was followed so that it was easy to add additional random network models. The random networks may either be created by randomizing real networks implementing the edge shuffling algorithm or by generating new random networks by using the available models (we used the Barabasi-Albert model [28] as the networks generated by using this model are scale-free i.e. they follow power-law degree distribution).

Network Topology Analysis
The topology of the network was analyzed using Net-workAnalyzer 3.3.1 [29] and CentiScape 1.2.1 [30] plugin fully implemented in Cytoscape 3.5.1. Various topological features calculated were: • Degree: Number of connections emerging from a node. It is one of the most elementary characteristics of a node. It is expressed as: Where, D i is the degree of node D and X ij is the adjacency matrix whose value is 1, if there exists an edge between the nodes i and j. If there is no edge between the two nodes then its value is zero. In our study, we have used this parameter to identify hubs and bottlenecks.
• Diameter [31] is the measure of the distance between the two most distant nodes. It is indicative of compactness of the network. A high value of diameter would suggest that the graph is not compact with respect to the nodes being considered. So, the value of diameter is a function of the two nodes that are chosen for distance calculation. Therefore, we see that a high value of diameter is not so much reliable as compared to a low value of diameter because such a value suggests that the graph is compact. • Density: It has been observed that biological networks are mostly sparse thus preserving robustness. The property used to calculate how dense or sparse a network is called network density. • Clustering Coefficient, (C) [31] is the measure of the cluster formation tendency of the graph. A cluster can be visualized as a subset of the vertices of the entire graph having a large number of edges of a node having 'm' neighbors is calculated by the formula: [31] of a node is the average connectivity of all the neighbors where the connectivity of a single node is given by the number of neighbors of that node. • Hub is a highly connected node and has a critical function in the network. Hubs are the nodes that have functional significance. • Hub Bottlenecks are nodes having high degree and betweenness centrality values and are significant in scalefree biological networks. Therefore, these nodes are termed as Hub-Bottlenecks (HBNs) [32]. These nodes are likely to be essential to the system and often coincide with high degree hub nodes [33]. • Closeness Centrality is the reciprocal of the summation of shortest distances between a vertex and every other vertex present in the network [34].
It is calculated as, is the shortest distance between node i and j.
• Stress Centrality: For a node 'a', the stress centrality is the number of shortest paths passing through 'a'. Higher the value of stress centrality the greater is the importance of a node in holding various communicating nodes together. It can thus help us identify connecting proteins. A node with a high-stress centrality value can potentially be a good target. • Topological Coefficient is a relative measure of the tendency of the nodes to share neighbors with other nodes. • Eccentricity is the minimum distance between a node, say b, and all other nodes. It is an index of the centrality of the node. In biological networks, nodes with a higher value of eccentricity are functionally accessible by other nodes of the network, and hence the changes in the concentration of the molecules are reflected readily. A protein having an eccentricity greater than the average eccentricity of the network has the capacity to influence other proteins.

Functional Annotation and Disease Gene Identification
KOBAS 3.0 [35] webserver was employed to ascertain the overexpressed genes and disease pathways in the network. In order to identify the statistically relevant hits, hypergeometric test was applied and p-value cut of ≤ 0.05 was kept as the threshold. KOBAS 3.0 is a web server for functional annotation of proteins or genes along with their functional enrichment. It integrates information from 9 gene set enrichment (GSE) methods. These methods can either be net-based or set-based. More than 4000 species and multiple pathways databases are supported by KOBAS 3.0. It supports both microarray and RNA-Seq data.
The results from KOBAS 3.0 helped in identifying the overexpressed pathways and genes. These pathways were superimposed on the integrated network. This helped in visualizing the genes involved in the pathogenesis of cancer both directly as well as indirectly.

Selection of the Novel Drug Targets from the Network
The approved drug targets in case of cancer epigenetics were searched and their interactions, as well as other topological characteristics in our network, were analyzed. The cytoHubba [36] plugin was utilized to rank the nodes and identify which nodes are capable of being classified as hubs. This plugin uses 11 topological methods among which we selected maximal clique centrality (MCC) owing to its better performance in terms of results and evaluation. The MCC of approved drug targets and rest of the network proteins were compared to identify potential drug targets.
For any node a, the MCC of a is defined as: S (a) is the collection of maximal cliques which contain v, and (|C|-1)! is the product of all positive integers less than |C|. If there is no edge between the neighbors of the node a, then MCC (a) is the same as the degree [37].

Cancer Epigenetics Network Structure
In order to develop and study HCC epigenetic pathway, a network consisting of all the DEG was developed. These DEGs were acquired from publicly available microarray study data of normal versus cancerous subjects. 250 DEGs were selected from each of the datasets i.e. GSE18081 [19], GSE37988 [20], GSE44970 [21], GSE54503 [22] and GSE57956 [23] expression profile using the cutoff of pvalue < 0.05 and |log FC| value >1.5. Furthermore, these DEGs were integrated with pathways involved in regulation of HCC. PPI network of DEGs and those searched from literature was constructed using STRING database [25]. Both the PPI networks were merged using the merge tool of Cytoscape 3.5.1 to develop the final epigenetics network consisting of 1423 nodes and 5282 edges (Fig. 3).

Network Validation
The network validation was conducted with the help of Network Randomizer plugin of Cytoscape 3.5.1. As compared to random network, cancer epigenetic networks have a higher value of average clustering coefficient depicting its modular nature. In our case, the average clustering coefficients for the random and epigenetic networks were 0.004 and 0.609 respectively. Besides, a higher value of average clustering coefficient also indicates redundancy and cohesiveness of the neighbors in the network. Random graphs have a small value of average clustering coefficient. Since proteins do not work in isolation, therefore, high cluster forming tendency of epigenetic network nodes indicates the biological significance of interactions between nodes. Secondly, the epigenetic network followed Power Law (Fig. 4) with a degree exponent 'γ' value of -0.778 and R 2 value of 0.817. Small 'γ' value denotes essentiality of highly connected nodes in the network and R 2 value closer to 1, indicating strong correlation between network nodes [38].

Topological Characteristics of the Network
Network characteristics as calculated from Network Analyzer 3.1.1 and Centiscape 2.2 are given in Table 1.
CytoHubba [36] plugin of Cytoscape 3.5.1, aided in identifying the central element of the cancer epigenetic network. Hubs are those proteins that interact with many partners and occupy the central region of a PPI network. HBNs are important linkers as they connect the sub-networks and have a high value of betweenness centrality. In our network, 25 and 13 nodes were identified as hubs and HBNs, respectively. Furthermore, the shortest path length distribution of the epigenetics network (Fig. 5A) also showed that the network is scale-free. This metric is widely used for analysis of diseaseassociated genes as functionally related genes tend to remain closer to each other. Similarly, the shared neighbor distribution (Fig. 5B) showed that the nodes in the network do not show the tendency of forming shared clusters rather every cluster is attached with each other by means of a single gene. These connecting genes are very important in the context of being considered as novel targets as these are crucial links between various pathways which eventually lead to epigenetic alterations [39].
Among these overrepresented pathways, 32 pathways have already been reported for their role in HCC epigenetics.
The novel ones identified are NOD like receptor signaling pathway, Fc epsilon signaling pathway and Hippo signaling pathway.

Identification of Potential Drug Targets for Cancer Epigenetics
The identification of potential drug targets was done on the basis of MCC. Therefore, the network nodes were ranked from the top 700 nodes which were selected based on their MCC score. From the list of 700 genes, we performed a rigorous literature survey to find out which of these genes have not yet been reported for their role in cancer epigenetics. We finally narrowed down our list to 10 genes based on their topological properties (Table 2, Fig. 6) which play a major role in HCC epigenetics and thus prove to be potential for epigenetic therapy. Furthermore, we also searched targets for investigational drugs at various stages of clinical trials used for cancer treatment and calculated their MCC values (Tables 3 and 4) [40]. The comparison between MCC values of approved and potential drug targets revealed that the novel targets cdk6, cdk4, cdkn2a, smad3, smad7, ccnd1, e2f1, sf3b1, ctnnb1 and tgfb1 have MCC score greater than the approved drug targets. Comparison of Tables 2 and 3 showed that our potential drug targets had features similar to that of approved drug targets. The value of eccentricity, which measures the easiness of a protein to influence the activities of several other proteins was greater than 1 for both approved and potential drug targets. Furthermore, both approved and potential drug targets had almost similar degree except cdkn2a that had a degree of 21. Similarly, the clustering coefficient and radiality of both approved and potential drug targets were close to 1. The radiality value is a measure of functional relevance of a protein for several other proteins but with the possibility of functionally irrelevant for few other proteins. The value of topological coefficient that describe the tendency of a node to have shared neighbors was close to 1 for approved drug targets, while our potential drug targets had much smaller values for topological coefficient.
The present clinical significance of the genes reported here are: • cdk6 activates cell proliferation. Overexpression of cdk6 is associated with resistance to hormone therapy in breast cancer. Upregulation of cdk6 has been observed in almost one-third of medulloblastoma cases. Changes in cdk6 can effect various hallmarks of cancer like induction of angiogenesis and evasion of growth suppressors etc.  • smad7 exerts an inhibitory effect on the EGF signaling pathway which is involved in breast cancer and ovarian cancer invasion and metastasis. A mutation located in smad7 gene is a cause of susceptibility to colorectal cancer (CRC) type 3. It is also overexpressed in pancreatic cancer.
• smad3 plays an important role in the regulation of genes involved in growth differentiation and death implying that an alteration in its activity or repressing its activity can lead to the formation or development of cancer. It has an important role in colorectal, breast and pancreatic cancer.
• ccnd1 codes for the protein cyclin D1. The overexpression of this protein is strongly correlated to ER+ breast cancer and deregulation of cyclin D1 is associated with hormone therapy resistance in breast cancer.
• e2f1encodes protein that play a crucial role in the controlling cell cycle and tumor suppressor protein. Elevated expression of E2F1 protein in breast cancer cell lines and head and neck carcinoma cell lines, and overexpression of E2F1 in invasive ductal breast cancer and non-small cell lung cancer. • sf3b1 gene mutations have been recurrently seen in cases of advanced chronic lymphocytic leukemia, myelodysplastic syndromes and breast cancer. sf3b1 is one of several genes involved in RNA splicing that has been identified as recurrently mutated in MDS and other malignancies. • ctnnb1 gene encodes for the protein β-catenin. This protein has dual function of managing transcription and intercellular adhesion. Mutations and overexpression of this protein has been found to be associated with many cancers like endometrial cancer, ovarian cancer and colorectal cancer.   • tgfb1 gene encodes for multifunctional peptide set TGFβ. It acts negative autocrine growth factor and controls the immune system response. Besides, it performs functions like cell growth, differentiation and apoptosis.

CONCLUSION
Graph theory has proven to be a promising tool in understanding various biological networks in the past. In the current work, this powerful tool implemented by means of Cytoscape 3.5.1 was analyzed to understand HCC epigenetic pathways. The network approach provides a deeper insight into the interactions taking place.
Epigenetic therapy is still in its infancy with only a few epigenetic drugs available. The novel targets reported here have the potential for opening new avenues in HCC therapy. Some of the reported targets are not directly involved in the pathways yet they play an indispensable role. One of the major reasons for interest in developing drugs based on epigenetics (epidrugs) is the reversible nature of epigenetic changes.
The major challenges epigenetic drug research is facing today and likely in the future prominently include interpreting the relevance of pharmacodynamics and somatic mutation level in the patient to hone the response to these agents.
Epigenetic alterations can lead to various other diseases like cardiovascular disorders, metabolic disorders, and neurological disorders. Since epigenetic changes are reversible, a better understanding of these mechanisms can be utilized to regulate the genetic switch and bring back the cell to the normal state.

ETHICS APPROVAL AND CONSENT TO PARTICI-PATE
Not applicable.

HUMAN AND ANIMAL RIGHTS
No Animals/Humans were used for studies that are the basis of this research.

CONSENT FOR PUBLICATION
Not applicable.

FUNDING
None.