Network properties of cancer prognostic genes and their gene sets in the human protein interactome

Background Identifying prognostic genes (PG) is crucial for estimating survival time and providing pinpoint treatments for patients with cancer. However, prognostic genes sets (PGS) reported in most existing research have low reproducibility and overlap ever between the same cancers or their subtypes. Their common characteristic as well as the molecular mechanism of action is still elusive. Methods Here, we obtained nine prognostic gene sets (including 1,439 prognostic genes) of different types of cancer from 23 high quality literatures, and systemically investigated eight network topological properties for PG and PGS compared with background and four other gene sets (cancer gene set CA, essential gene set ES, housekeeping gene set HK, and metastasis-angiogenesis gene set MA) based on the HPRD and String networks. Results The results showed that PG did not occupy key positions in the human protein interactome network, and were more similar to ES rather than CA. Also, PGS had significantly small intraset distance (IAD) and interset distance (IED) in comparison with random sets. Further, we also found that PGS tended to have be distributed within network modules rather than between modules, the functional intersection of the modules enriched with PGS was closely related to cancer. Conclusions Our research reveals the common properties of cancer PG and PGS in the human protein interactome network, and can help us understand and discover cancer prognostic biomarkers.


Abstract
Background Identifying prognostic genes (PG) is crucial for estimating survival time and providing pinpoint treatments for patients with cancer. However, prognostic genes sets (PGS) reported in most existing research have low reproducibility and overlap ever between the same cancers or their subtypes. Their common characteristic as well as the molecular mechanism of action is still elusive. Methods Here, we obtained nine prognostic gene sets (including 1,439 prognostic genes) of different types of cancer from 23 high quality literatures, and systemically investigated eight network topological properties for PG and PGS compared with background and four other gene sets (cancer gene set CA, essential gene set ES, housekeeping gene set HK, and metastasis-angiogenesis gene set MA) based on the HPRD and String networks. Results The results showed that PG did not occupy key positions in the human protein interactome network, and were more similar to ES rather than CA. Also, PGS had significantly small intraset distance (IAD) and interset distance (IED) in comparison with random sets. Further, we also found that PGS tended to have be distributed within network modules rather than between modules, the functional intersection of the modules enriched with PGS was closely related to cancer. Conclusions Our research reveals the common properties of cancer PG and PGS in the human protein interactome network, and can help us understand and discover cancer prognostic biomarkers.

Background
Prognostic genes (PG) have many crucial clinical applications, such as accurate predictions of cancer types (or subtypes), stages and their survival time for cancer patients. In particular, precise targeted treatments and surveillance strategies could be implemented when patients have been classified into different risk groups by means of 3 application of PG [1]. In the past 20 years, tremendous efforts have been making to investigate PG, and a large amount of prognostic biomarkers have been identified for different types of cancer [2][3][4][5][6][7][8][9][10], Several PG have been playing critical roles in the prognosis of certain types of cancer such as ER and HER2 for breast cancer [11].
Biological network provides an invaluable research platform to trace genetic phenomena and disease mechanisms on a system level [12,13]. Network topology analysis helps to discover groups of nodes with special network characteristics in biological networks, as well as associations between groups (e.g., plant immunity[14] and human disease [15,16]). Among them, study of cancer genes network has showed that they tend to have higher degree and betweenness compared with essential genes [17,18]. Systematic studies of the PG's network properties could help identify pan-cancer PG and unveil the mechanisms of cancer. However, previous studies of PG were scattered and mostly focus on one specific cancer type or subtype. Cumulative evidence also showed that even for the same cancer type, PGS obtained by different researchers had very small overlap and questionable reproducibility [19]. Few PG studies were carried out involving multiple cancer types (e.g., [20][21][22] ). Also, they either didn't pay attention to the network properties of PG, or were only involved in the gene co-expression network characteristics of PG in a few cancers. So, it is still not well known to their common properties of the human protein interactome.
In this study, we obtained nine PGS (including 1,439 PG) of different cancer types when we manually selected 23 high quality literatures that involve the identification of prognostic genes. Based on the two protein interactome networks (HPRD and String) and four other gene sets for comparison (cancer gene set: CA, essential gene set: ES, housekeeping gene set: HK, and metastasis-angiogenesis gene set: MA), we investigated their eight network properties. Our study showed that although PG did not possess a high degree compared 4 with cancer genes, PGS had tighter network connections and closer inter-genesets distances than background, and the network modules they were in had many common functions that were closely related to cancer. These findings could bring insight into identification of prognostic biomarkers and understandings of their mechanisms.

Results
Overview of prognostic genes. To systematic study of prognostic genes, we obtained 25 small prognostic gene sets (ranging from 3 to 330) from 23 high quality literatures (Table 1). Similar to previous study [19], these genes had very small overlap and network connections. Only 14 genes were repeatedly mentioned 3 times in these small gene sets (see Supplementary Table 1 Table 2). To investigate their network properties, we employed two protein-protein interaction (PPI) networks, HPRD and String, which both exhibited power-law node-degree distributions ( Figure 1A and B) [23]. Figure 1C shows that cancer prognostic genes are discretely distributed in the HPRD network. Only three out of the 14 genes which appeared three times above had directly connected edges in the HPRD network.  Figure 2. Like ES, degree and betweenness of PG was lower than the background, while CA and MA were obviously higher than the background in two PPI networks in Figure 2A-D. However, in Figure 2E-H, Closeness of PG and four other gene sets were considerably higher than the background, while Eigenvector of PG was different from CA and MA, and its values were always lower than the background in the HPRD and String networks. Eigenvector of CA and MA, as well as degree and betweenness of HK, showed inconsistency in the both networks, probably due to the String network consists of more notes and edges [24].
Overall, the results clearly showed that: (1) CA had very similar centralities to MA and both gene sets had significantly higher centralities than other three gene sets including PG (except eigenvector in the String network, FDR-adjusted p-values of t tests were much smaller than 0.001 in all other cases). This illustrated that PG was significantly different from cancer-related genes in terms of four network centralities; (2) Except Closeness, the other three centralities of PG were below the average of the whole network. This showed that, unlike cancer genes, prognostic genes did not occupy key positions in the network [18,25]; (3) The four centralities of PG were not significantly different from those of ES (FDR-adjusted p-values of t test were greater than 0.1 for all four centralities).

Four network measures of prognostic gene sets. Most of cancer prognostic
biomarkers often act as functional units in a gene set. Therefore, it was necessary to examine the network topological properties of gene sets. To this end, we first calculated clustering coefficient (CC) for nine PGS, four other gene sets and random gene sets. CC measures the tendency that the nodes in a graph cluster together. Larger CC values indicate that the nodes are more likely to form clusters in a network [26]. Figure 3A and B shows their distributions of CC. In the HPRD network, nine PGS had slightly larger CC than 6 the random gene sets (p-value of KS test was not significant). CA and MA also had larger CC than HK and ES. For the String network with higher density, on the one hand, nine PGS had significantly smaller CC than the random gene sets (KS test, p-value < 0.05), which showed that genes within the nine PGS were more sparsely connected compared to random gene set in the network. On the other hand, HK had significantly larger CC than all the other gene sets (p-value of permutation test smaller than 0.001). This was probably due to the fact that edges were more likely to be formed between HK in the String network since HK had consistent expression patterns [27]. And it can also be clearly demonstrated by comparing the degree of HK in the two networks (Figure 2A and B).
Through the investigation of CC, we failed to obtain the significant common properties of the PGS in the network. So, we proposed three other measures, intraset distance (IAD), interset distance (IED) and genset-distribution in modules (GDM), to examine the network properties of gene sets in the network. IAD and IED were used to portray the network distance within a gene set and between two gene sets, respectively (see the methods section for more details). Their calculations were based on shortest path (SP), which can reflect the ability of network information transfer [28]. Figure 3C and D shows that IAD of nine PGS are significantly smaller than the random gene sets in both networks, indicating that there is a more compact network structure within PGS. In the four other gene sets, CA and MA had obviously smaller IAD than PGS compared to HK and ES, and considering two networks together, the ES was not the closest one to PGS in the IAD distribution.
Similarly, we found that IED between PG themselves were significantly smaller than those between PGS and the random gene sets ( Figure 3E and F). The result indicates that the nine PGS are not spatially loose, but rather closely connected. Simultaneously, we also found that IED between PGS and the four other gene sets were significantly smaller than those of PGS themselves. One possible reason was that they were derived from different 7 cancer types. Among them, we found that IED between PG and CA or MA was smaller than the other two gene sets, which may suggest that PGS is more closely related to cancer ( Figure S1A and B). In addition, we can easily see that whether it was IAD and IED, the distances in the String network were smaller compared to the HPRD network.
Next, we used genset-distribution in modules (GDM) to investigate the distribution of PGS within and between modules of network. Figure 3G and H shows that GDM of nine PGS are significantly larger than random in both networks, demonstrating that they are more likely to be distributed within the module. We also found that GDM of MA were the largest of the four other gene sets, and a change in the relative position of CA and HK in the two networks. This may be due to the different modules that were derived from different networks, and the complexity of gene sets themselves.
Functional analysis of prognostic gene sets. We performed functional enrichment analysis for nine PGS based on GO terms and KEGG pathways database using Fisher test.
However, more than half of the gene sets were not enriched with any significant functional terms. Genes with same or similar functions are more inclined to be in the same module of a network [29]. We then examined functions of network modules with two or more PGS.
Interestingly, when we compared these functions of the modules from different networks, we found that the intersections of their functions were mostly related to cancer. Figure 4 shows the intersection of the function of module #4 of the String network and module #5 and #7 of the HPRD network. Most of functional terms could be attributed to hallmarks of cancer [30]. They included "Extracellular matrix organization", "Leukocyte migration", "Collagen metabolic process", "Transforming growth factor beta receptor signaling pathway", etc. In particular, among them, "extracellular matrix organization"was the most significant GO term. Researchers have found that its remodeling directly affects tumor growth, development, and progression [31]. And "transforming growth factor beta" (TGF-β), the main pathway of another functional terms, have been evaluated as prognostic or predictive markers for cancer patients [32]( the lower half of Figure 4).

Discussion
A few nodes, the hubs, have a higher connectivity coexistence with most rarely connected nodes in a scale-free network, and they are highly influential in keeping the whole network together [33]. Hub genes were found in human cancer genes and essential genes of yeast and worms in PPI networks [17,23]. In the study of network centralities, we also found that the degree of CA and their other three centralities were significantly larger than PG and human ES. Among them, PG and human ES were very similar (their distribution difference was not significant by t-test) in most cases, and were smaller than the background mean state, although ES was reported to have significantly more important network topology than "unnecessary genes" [34]. However, low-degree features did not affect functional genes to play an important role. For example, the metabolites with low degree were involved in essential reactions in the metabolic networks of Escherichia coli Although the number of genes in prognostic biomarkers had a decreasing trend as a whole [39,40], most of the previous prognostic biomarkers were often in the form of a union set of dozens or hundreds of genes, or even a network module [6,[41][42][43]. This would considerably weaken the importance of individual genes, which may be one of the reasons why the four network centralities of the single PG were not high. However, by focusing on gene sets, our study helped to make up for this deficiency. Three network properties of gene sets presented in this paper had obtained results that were significantly different from the random background. They were more conducive to further understanding of prognostic biomarkers and their mechanisms of action. Interestingly, Yang et. al. [22] also found that prognostic genes did not occupy hub positions and were more likely to appear within network modules when studying the topological properties of prognostic genes in co-expression networks based on four types of cancer. Despite using the different measurement method, Zhang and Horvath [44] also drew similar conclusions that prognostic genes for cancer survival was highly correlated with their intramodular connectivity. The network structure of PPI may be less susceptible to environmental conditions than gene co-expression network [16]. This implies that the topological properties of prognostic genes in protein interactome networks can also be shared in gene co-expression network.
Discovery and clinical applications of prognostic genes have been going on for several decades, but the secrets of prognostic genes are yet to be unveiled. Non-overlapping and non-reproducibility of research findings on prognostic genes can be due to many factors, such as types of cancer, microenvironment, patient cohorts, methodologies and technological platforms [20]. Although there are enormous difficulties to the systemic study of PG, further studies will be still essential in view of the indispensable roles of prognostic genes in cancer diagnosis and treatments as well as the exploration of cancer mechanisms. The above unfavorable situation may be improved by studies of emerging new molecular markers, such as, prognostic miRNAs, prognostic lncRNAs, prognostic circRNAs, and their combinations [45][46][47].

Conclusions
In summary, we systematically studied the network properties of PG and their PGS for the first time based on two protein interactome networks and eight network properties. We found that prognostic genes had significantly different network properties from CA and were similar to ES in the four network centralities. For intra-module and inter-module distances, PGS was significant smaller relative to random gene sets. And they were more easily enriched inside the modules, which were found to enrich the functions related to cancer development and progression. These characteristic will be valuable for future understanding and identification of prognostic genes. However, several disadvantages still existed in this study, including not considering emerging markers, lacking optimized combinations of gene sets, fewer gene sets and cancer types, etc. Including more datasets and developing new computational strategies could lead to more significant results in the future research.

Methods
Prognostic genes and other four gene sets. To obtain reliable prognostic genes from vast amounts of existing literature, we set two criteria for a publication to be included in our study: (1) the sample size of patients with cancer in the study was larger than 100; (2) the publication was published on a high-impact biomedical journal (impact factor > 8.0) after year 2000 and the research findings are explicit and accessible. Given these criteria, we manually selected 23 publications and each publication reported 3 to 330 prognostic genes. For convenience we also merged these genes into nine gene sets according to cancer types and the number of prognostic genes reported in each study. Each gene set consisted of 100 to 200 prognostic genes. Details of the selected publications and the corresponding prognostic genes can be found in Table 1  Clustering coefficient (CC) and shortest path (SP) were calculated using previous definitions [12]. Specifically, if two nodes were not connected in the network, their SP was set to be the maximum of SP between all connected nodes in the network.
We proposed two measures to further quantify the network properties of gene sets. First, we used SP to define the distance between two nodes in the network. Then, we defined intraset distance (IAD) and interset distance (IED) based on the distance. IAD was derived from the definition of the average shortest path of complex network [53]. These two measures were used to quantify the distance (or compactness) within a gene set and the distance between two gene sets, respectively. For gene set S with N genes, the IAD was defined as follows: 12 Here, E intra is the total number of edges within network modules and E inter is the total number of edges between network modules. Figure 5 demonstrates how GDM is calculated.    Schematic diagram of calculating GDM of gene sets in the network. The formula for GDM and its calculation process were provided for the given example in the chart below.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.