BACKGROUND
researchers often use pathway analysis  <cit>  to reveal novel insights into their gene lists obtained from high-throughput experiments.
one simple approach is evaluating the number of genes from a differentially expressed gene list found in a particular pathway.
advanced methods calculate pathway-level statistics.
in some studies, a pathway is also considered to be a gene set, a group of genes sharing common biological functions.
for example, gene set enrichment analysis  calculates gene set-level statistics  <cit> .
using gene set-level statistics has additional advantages because significant analysis at the single gene level suffers from a limited number of samples and noise  <cit> .
gene set based methods have also been developed to investigate phenotypic changes at the pathway level  <cit> .
by using pathway topology, researchers can obtain a better ranking of pathways/gene sets from their gene lists  <cit> .
several approaches have been proposed for constructing networks of pathways/gene sets to study complex molecular networks.
yong, et al.
constructed a global pathway crosstalk network and linkage network for yeast  <cit> .
dikla, et al.
proposed a method for gleaning patterns of interactions among biological processes by analyzing protein-protein interactions, transcriptional co-expressions, and genetic interactions  <cit> .
gene set networks can be also applied to understanding diseases.
liu, et al.
proposed an approach to detect the crosstalk among alzheimer's disease  related pathways and the dysfunctions in the six brain regions of ad patients  <cit> .
for databases of pathway/gene set network, sudhir, et al.
developed human pathway database  to enable the study of human pathway networks  <cit> .
huang, et al.
developed pathway and gene-set enrichment database  to further study human gene set networks  <cit> .
another study proposed an approach to examine interactions between pathways in mice by integrating different types of data  <cit> .
recently, jignesh, et al.
proposed methods to construct multi-edge gene set networks  to reveal insights into global relationships among biological themes or gene sets  <cit> .
a multi-edge gsn consists of three types of edges: co-membership edge , linkage edge , and co-enrichment edge .
to our knowledge, none of the existing gene set or pathway networks contain high-level directionality information, which is essential for building models to understand the mechanistic linkages among many gene sets or pathways that are enriched in a biological condition.
in this paper, we present a method of constructing a directed, regulatory gene set network called regulatory gene set network .
a directed edge in an r-gsn represents a regulatory relationship from one gene set to another.
compared to a linkage gene set network , which is constructed from protein interactions discarding directionality of the interactions, an r-gsn provides higher resolution of knowledge.
our hypothesis is that r-gsns can reveal novel gene set relationships and provide complementary knowledge to the existing types of gsns, such as m-gsns, l-gsns, and e-gsns.
in an r-gsn, a pair of gene sets are connected if a significant number of gene regulations exist between the unique genes of the gene sets.
gene set and gene regulation data were collected from multiple sources.
to evaluate r-gsns, m-gsns are constructed separately as baselines for each collection including kegg, reactome, and three types of gene ontology  terms.
we exploited the directionality information provided by r-gsns to search for significant gene sets in the network.
we also compared our r-gsns with e-gsns obtained from jignesh et al  <cit> .
we chose alzheimer's disease  as an example to study disease-specific r-gsns.
in order to construct an ad specific r-gsn, gene sets from different collections were first combined into one single collection and then ad related gene sets were used to construct the ad specific r-gsn.
methods
in a gene set network , a node represents a gene set and an edge represents a relationship between two gene sets.
a gsn helps explain biological complexity by revealing high level relationships among biological processes.
we constructed a new type of gsn, a regulatory gene set network , using public available gene regulation data.
another type of gsn we constructed is the co-membership gene set network , which can be a baseline network as it is constructed from annotated gene sets and thus provides experimental validation.
finally, we used hypergeometric distribution to calculate significance values for each edge in both r-gsns and m-gsns.
data sources
the two major types of data used to construct r-gsns are gene set data and gene regulation data.
gene set data provides information for nodes in a gene set network as well as information for constructing co-membership gene set network.
five collections of gene sets were collected.
they included kegg, reactome, go biological process, go cell component, and go molecular function from kegg  <cit>  and msigdb  <cit> .
the total number of gene sets from the five gene set collections is  <dig>  and the total number of genes is  <dig> .
gene regulation data was used for constructing r-gsns.
human gene regulations from string  <cit> , transfac  <cit> , tred  <cit> , and spike  <cit>  were collected.
different types of id of genes obtained from different data sources were mapped to ncbi official gene symbols.
in order to select only the high quality human gene regulation data, different criteria were used to filter gene regulations in different data sources.
for string, gene regulations with score >=  <dig> were collected; for transfac, binding site quality <=  <dig> was the criterion; and for tred, only gene regulations not obtained from computational predicted method were collected.
all gene regulations provided by spike were collected since they are from pathway.
the total count of unique gene regulations after combining data from the four data sources is  <dig> .
construction of gene set networks
an m-gsn was constructed separately for each gene set collection.
the hypergeometric distribution was used to calculate the significance value for each co-membership edge.
<dig>  count the number of genes inside each of the two gene sets, gs <dig> and gs <dig>  and the number of shared genes between gs <dig> and gs <dig>
<dig>  calculate the p-value by using the hypergeometric distribution:
p-value=kkn-kn-k/nn
where n is the total number of genes; n is the number of genes in gs1; k is the number of genes in gs2; and k is the number of shared genes.
<dig>  adjust the p-value for multiple hypotheses to control the false discovery rate by using the benjamini-hochberg procedure with p-value ≤  <dig> .
<dig>  connect a pair of gene sets with an edge if the edge is rejected by the benjamini-hochberg procedure and thus considered as a significant edge.
an r-gsn was constructed as follows:
<dig>  count the number of gene regulations where genes in gs <dig> regulate genes outside gs <dig> and the number of gene regulations where genes outside gs <dig> regulate genes in gs <dig>
<dig>  remove shared genes from gs <dig> and gs <dig> and count the number of gene regulations where the remaining genes in gs <dig> regulate the remaining genes in gs <dig>
<dig>  calculate the p-value by using the hypergeometric distribution using formula , where n is the total number of gene regulations; n is the number of gene regulations where genes in gs <dig> regulate genes outside gs2; k is the number of gene regulations where genes in gs <dig> are regulated from genes outside gs2; and k is the number of gene regulations from genes in gs <dig> to gs <dig>
<dig>  adjust the p-value for multiple hypotheses to control the false discovery rate by using the benjamini-hochberg procedure with p-value ≤  <dig> .
<dig>  connect a pair of gene sets with a directed edge pointing from gs <dig> to gs <dig> if the edge is rejected by the benjamini-hochberg procedure and thus is significant.
r-gsns were constructed separately for each gene set collection.
note that it is possible that r-gsns contain a loop when both incoming and outgoing edges of a pair of gene sets are significant.
r-gsns are directed networks, while m-gsns are undirected networks.
before comparing r-gsns with m-gsn, the r-gsn was converted to an undirected network by discarding the directions of edges and removing loops.
a disease specific gene set network
a disease specific gene set network was constructed in order to show how r-gsns can help researchers explain disease complexity.
alzheimer's disease  was chosen to be a case study.
an ad gene list containing  <dig> genes from alzgene database  <cit>  was first collected.
the five gene set collections were combined into a single global collection without removing and changing the original gene sets.
from the list, ad related gene sets were selected from the global gene set collection.
then the ad specific gene set network was constructed.
in the network, each node represents an ad related gene set.
in order to find ad related gene set, the number of genes in each gene set found in the ad gene list was counted.
each ad gene list was treated as a new gene set.
the same method used in constructing m-gsns was used to calculate the p-value for each gene set.
there were  <dig> gene sets which shared significantly high numbers of genes with the ad gene list.
these ad related gene sets were collected and used to construct the ad specific r-gsns.
network analysis
three types of centrality  <cit>  were calculated for each gene set in both r-gsns and m-gsns, degree centrality, betweeness centrality, and closeness centrality.
igraph software package for r  <cit>  was used in order to compute all network values.
degree centrality represents the number of edges upon a node.
a gene set with a high degree centrality is likely to be an important gene set because it acts like a hub in the network.
it can also be used for comparing two different types of networks.
to conduct such comparison, we calculated pearson's correlation coefficient of nodes' degree centrality between an m-gsn and an r-gsn.
betweenness centrality is defined as the number of times a node acts as a bridge along the shortest path between two other nodes.
in a gene set network, a gene set with a high betweenness value is likely to be a part of several biological critical paths.
closeness centrality is defined as the inverse of the average length of the shortest paths between a node and all other nodes in a network.
a node with a high closeness value is more central.
in a gene set network, if a gene set with a high closeness value is disturbed, it is likely that a high number of gene sets will be affected by the gene set.
RESULTS
interpreting co-membership gene set network and regulatory gene set network
when a pair of gene sets in an m-gsn are connected, the interpretation depends on gene set data types.
for pathway gene sets, an m-gsn represents pathway crosstalk; whereas for go gene sets, an m-gsn represents protein moonlighting or gene sharing.
when a pair of gene sets in an r-gsn are connected, they are connected by a directed edge.
a directed edge in an r-gsn presents a possibility of one gene set regulating another.
according to table  <dig>  the numbers of unique genes from all gene set collections are not much different, where reactome and go biological process contain more genes than the others.
the total count of unique genes from all collections is  <dig> .
we calculated the proportion of regulatory edges to nodes and normalized the proportion value by using number of unique genes.
we used rigi/g/si, where ri is the number of gene is sets in collection i; gi is the number of unique genes in collection i; g is the total number of unique genes ; and si is the number of gene sets in collection i, to calculate the normalized proportion.
in the r-gsn, go biological process has the highest normalized proportion of edges to nodes  among the five gene set collections.
these results indicated that pairs of biological process gene sets are more likely to have a regulatory relationship.
in addition, according to the distribution of gene set sizes , the five gene set collections have similar distributions.
the distributions show that the five collections contain more small gene sets  than large gene sets .
afor regulatory edges, a pair of gene set can have a loop if both incoming and outgoing edges are significant.
the number of regulatory edges is always greater than or equal to the number of regulatory relationships because one pair of gene sets either have or do not have regulatory relationship.
comparing the kegg regulatory gene set network and the kegg co-membership gene set network
in the m-gsn of kegg, an edge between pathways can be considered as a pathway crosstalk  <cit> .
therefore, m-gsns can be baselines to compare with r-gsns.
a regulatory relationship between pathways can provide complementary knowledge such as dysfunction of one pathway affecting function of other pathways.
to investigate this knowledge, r-gsns and m-gsns of the kegg pathways were compared.
the percentage of shared edges between m-gsns and r-gsns is the highest for kegg pathway  =  <dig> %).
the low percentage of shared edges indicates that the r-gsn provides complementary information to the m-gsn.
in addition, it is important to note that r-gsns are constructed from gene regulations which were collected from high coverage data sources.
therefore, it is unlikely that r-gsns depend on the quantity and the quality of experimental data.
considering both the r-gsn and m-gsn of kegg pathway gene sets, the most significant edge of the kegg r-gsn is a regulatory relationship from "cell cycle" to "cytokine-cytokine receptor interaction" with p-value =  <dig> e- <dig> ; whereas the co-membership edge between "cell cycle" to "cytokine-cytokine receptor interaction" has relatively low significance value .
in addition, only  <dig> of the top  <dig> most significant regulatory edges are found in m-gsns.
these findings suggest that the r-gsn reveals additional knowledge to the m-gsn.
the regulatory relationship is gene set  <dig> regulates gene set  <dig>
for the kegg r-gsn,  <dig> of the  <dig> most significant regulatory edges are from the "cell cycle" gene set to other  <dig> kegg pathway gene sets .
this suggests that changing of "cell cycle" pathway is likely to affect other pathways.
this finding is corresponding to the fact that a cell cycle is a complex series of phenomena by which cellular material is duplicated and divided.
therefore if a cell cycle pathway does not function appropriately, other pathways such as pathways in cancer can be affected.
for the m-gsn, the kegg pathway gene sets of alzheimer's, parkinson's, and huntington's diseases have significant co-membership edges linking them together .
the three co-membership edges connecting the three neurodegenerative diseases are among the top  <dig> most significant co-membership edges suggesting that the three neurodegenerative diseases are highly related.
in addition,  <dig> of the top  <dig> co-membership edges connect cancer related pathway gene sets.
these  <dig> edges connect the "pathways in cancer" gene set with  <dig> gene sets of cancers including "small cell lung cancer", "pancreatic cancer", "melanoma", "colorectal cancer", and prostate cancer".
in addition, degree centrality  of each node in both networks was calculated .
in the kegg m-gsn, the "pathways in cancer" gene set has the highest value of dc .
in the kegg r-gsn, the "cell cycle" gene set has the highest dc  and the highest out-degree centrality .
the gene set which has the highest in-degree centrality  is "pathways in cancer".
using directionality information to calculate in-degree and out-degree centrality of the kegg r-gsn, we can find the sink, "pathways in cancer", and the source, "cell cycle", gene sets.
after the dc value for each gene set was calculated, the correlation between dc of the kegg m-gsn and dc of the kegg r-gsn was calculated.
the correlation coefficient is  <dig>  and r-squared value is  <dig>  .
this result suggests that the gene set which is important in the m-gsn is likely to be important in the r-gsn and vice versa.
we found three interesting outliers.
pathway number  <dig>  which has regulatory dc =  <dig>  and co-membership dc =  <dig> , is "cell cycle" suggesting that the cell cycle pathway does not tend to share genes, but to have regulatory relationship with other pathways.
pathway number  <dig>  which has regulatory dc =  <dig>  and co-membership dc =  <dig> , is "pathways in cancer" suggesting that pathways in cancer shared high number of genes with several pathways and their genes also regulate the unique genes of other pathways.
pathway number  <dig>  which has regulatory dc =  <dig>  and co-membership dc =  <dig> , is "maturity onset diabetes of the young".
this pathway shares  <dig> genes with only "type ii diabetes mellitus pathway".
while the correlation between dc of the kegg m-gsn and dc of the kegg r-gsn is as high as  <dig> , the topologies of the networks are different .
this suggests that the two networks can be used to explain different phenomenon in biological systems.
constructing an exclusive r-gsn
we constructed a kegg exclusive r-gsn  that contains only exclusive edges in the r-gsn .
the correlation of dc between the kegg exclusive r-gsn and the kegg r-gsn is  <dig> .
the correlation of dc between the kegg exclusive r-gsn and the kegg m-gsn is  <dig> , while the correlation of dc between the kegg non-exclusive r-gsn and the kegg m-gsn is  <dig> .
these suggest that the exclusive r-gsn can reveal important gene sets that are not likely to be revealed by m-gsn.
the directionality information from a regulatory gene set network revealed "sink" and "source" gene sets in addition to "hub" gene sets, which can be revealed by constructing.
table  <dig> shows the top  <dig> highest out-degree centrality gene sets  and top  <dig> highest in-degree centrality gene sets  of kegg exclusive r-gsn.
"cell cycle", "p <dig> signaling pathway", and "tgf-beta signalling pathway" are among the top  <dig> highest out-degree centrality which can be the sources.
the top  <dig> highest in-degree centrality gene sets, "cell cycle", "hematopoietic cell lineage", and "cytokine-cytokine receptor interaction", are the sinks.
comparison of the kegg co-enrichment network and the kegg regulatory network
in the previous analysis, we used the kegg m-gsn as a baseline and compared it with the kegg r-gsn.
several relationships between gene sets were found exclusively in the r-gsn.
in order to validate both the kegg m-gsn and the kegg r-gsn, the two gsns were compared with the kegg co-enrichment network  obtained from the study of jignesh, et al.
<cit> .
to construct an e-gsn, they integrate experimental gene lists and link two gene sets if the unique genes of the two gene sets are consistently enriched together across many experimentally derived gene lists.
therefore, edges found in the kegg r-gsn should also be found in the kegg e-gsn.
the total number of edges found in both the e-gsn and the r-gsn is  <dig> , which accounts for  <dig> % of the total number of edges in the e-gsn.
the total number of edges found in both the e-gsn and m-gsn is  <dig>  which is equal to  <dig> % of the total number of edges in the e-gsn.
in order to calculate the significance value of the number of shared edges, the kegg e-gsn was compared with random networks.
we randomly generated  <dig>  networks using all the  <dig> gene sets from kegg.
to calculate the significance value of the number of shared edges between the kegg e-gsn and the kegg m-gsn, each of the  <dig>  random networks contains  <dig>  edges, which are equal to the number of edges found in the kegg m-gsn.
then fisher's exact test was used for calculating the p-value for the number of shared edges.
the p-value is <  <dig> e- <dig> .
to calculate the significance value of the number of shared edges between the kegg e-gsn and the converted kegg r-gsn, each of the  <dig>  random networks contains  <dig>  edges.
then fisher's exact test was used for calculating the p-value for the number of shared edges.
the p-value is <  <dig> e- <dig> .
the number of shared edges between the kegg e-gsn and the kegg r-gsn is significantly high.
this is corresponding to the fact that a pair of gene sets with strong regulatory relationship should be connected with a co-enrichment edge.
the number of shared edges between the kegg e-gsn and the kegg m-gsn is also significantly high.
this is also corresponding to the fact that a pair of gene sets with a high number of shared genes should be connected with a co-enrichment edge.
a disease specific regulatory gene set network
the number of shared genes between the ad gene list  and each gene set in the global gene set collection was counted.
out of the  <dig>  gene sets,  <dig> have significant number of shared genes.
among the  <dig> ad gene sets,  <dig> are from kegg;  <dig> are from reactome;  <dig> are from go molecular function;  <dig> are from go biological process; and  <dig> are from go cellular component.
figure  <dig> shows the degree distribution of ad-specific r-gsn constructed based on these  <dig> ad gene sets.
in addition, table  <dig> shows the top  <dig> ad related gene sets arranged by the significant values.
"genes involved in lipid digestion, mobilization, and transport" gene set has the highest significant value,  <dig> e- <dig>  the significant value was calculated by using formula  where there are  <dig> genes and  <dig> genes were found in the ad gene list.
the ad related r-gsn of the top  <dig> gene sets are presented in figure  <dig>
we then investigated the top  <dig> degree centrality gene sets.
table  <dig> shows that "signal transduction" from go biological process has the highest value of dc suggesting that the signal transduction process is very important in ad.
searching on pubmed found more than  <dig>  publications discussing the relationship between signal transduction abnormality and alzheimer's disease.
considering the directionality of the ad specific r-gsn, "signal transduction" has the highest in-dc .
gene sets which have the highest out-dc  are "cellular protein metabolic process" and "cellular macromolecule metabolic process" .
these results show that directionality information from the r-gsn enables us to identify, in ad context, a sink gene set, "signal transduction", and source gene sets, "cellular protein metabolic process" and "cellular macromolecule metabolic process".
furthermore, the closeness and betweenness of each gene set in the ad specific r-gsn were computed.
note that the directionality of each edge in the network was considered when we calculated the shortest path for the network in order to calculate the closeness and the betweenness.
for the closeness centrality, "signal transduction" from the go biological process, "pathways in cancer" from kegg, and "protein metabolic process" from the go biological process have the highest closeness values .
these results suggest that inappropriate functions of the three gene sets are likely to affect high number of gene sets in ad.
for the betweenness centrality, "pathways in cancer" from kegg, "system development" from go biological process, and "leishmania infection" from kegg have the highest betweenness values .
these results suggest that the three gene sets are likely to be on the critical path of biological functioning for ad patients.
CONCLUSIONS
co-membership gene set networks  and regulatory gene set networks  for the five different gene set collections were constructed, compared, and studied for their biological relevance in this study.
the results show that new r-gsns can provide complementary biological information to conventional m-gsns.
the results also show that while the correlation between the degree centrality of the kegg m-gsn and the degree centrality of the kegg r-gsn is relatively high, the topologies of the networks are totally different.
this suggests that the two networks can be employed to explain different phenomenon in biological systems.
validating kegg m-gsns and kegg r-gsns separately against the kegg co-enrichment networks   <cit>  reveals that the numbers of shared edges are significantly high.
the results suggest that the use of e-gsn to validate m-gsn and r-gsn could be a better way to describe experiments.
moreover, the r-gsn specific to alzheimer's disease shows that the network can be used as a good high-level mechanistic model complementary to a more difficult-to-derive gene regulatory network model towards a systematic understanding of the disease mechanism.
in this study, we describe methods for constructing both global and disease specific r-gsns.
they enable future research on systems biology and translational bioinformatics.
since the underlying gene regulation data are collected from high quality and high coverage data sources, the directed edges in the r-gsn do not tend to depend on the number and the quality of experimental data.
moreover, directionality information from the r-gsn enables the finding of source gene sets and sink gene sets, which might be important for drug discovery or drug repositioning.
tissue-specific and disease-specific gene regulations were not used in this study.
therefore, the gsns obtained might be generally applicable but should be carefully reviewed and curated when they are to be used in specific disease contexts where differential gene regulations exist extensively.
however, this limitation can be alleviated by constructing only condition-specific gsns instead of global gsns in the future.
the framework in this study can also be later extended by collecting higher resolution such as tissue specific and disease specific gene regulation data and gene set data.
separately, we have been developing a comprehensive database populated with new gene sets and all relationships identified or integrated between genes and gene sets into an online resource called pager  <cit>  to put the method developed in this work into a web application for biological users to explore.
list of abbreviations used
gsn: gene set network; m-gsn: co-membership gene set network; e-gsn: co-enrichment gene set network; l-gsn: linkage gene set network; r-gsn: regulatory gene set network; dc: degree centrality; m: co-membership edge; e: co-enrichment edge; r: regulatory edge; l: linkage edge.
competing interests
the authors declare that they have no competing interests.
authors' contributions
jake y chen conceived of this work, guided the research team by providing ideas and feedback along the way, and revised the manuscript.
chayaporn collected the data, constructed gene set networks, analyzed the networks, and wrote the manuscript.
liugen assisted network visualization and revised the manuscript.
supplementary material
additional file 1
distribution of the sizes of gene sets in each collection.
this file contains five histograms presenting distribution of the sizes of gene sets in each of the five gene set collections.
click here for file
additional file 2
ad gene list.
this file contains a list of ad related gens.
click here for file
additional file 3
centrality values of each gene set in the ad specific gene set network.
this file contains centrality values, including degree centrality, closeness centrality, and betweenness centrality, of  <dig> ad related gene sets in the ad specific r-gsn.
click here for file
acknowledgements
we thank dr. xiaogang wu for his generous help in conceptualizing the work, providing occasional comments, and proofreading the manuscripts.
we thank the support from indiana university information technology services for their partial support of this work through high-performance computing hardware and software.
we thank indiana center for systems biology and personalized medicine for providing research environments that led to the timely completion of this work.
declarations
publication costs for this article were funded by the corresponding author's affiliated institutions, indiana university school of informatics and computing in the united states and wenzhou medical university in china.
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2015: selected articles from the fourth ieee international conference on computational advances in bio and medical sciences : genomics.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/16/s <dig>

 researchers often use pathway analysis to reveal novel insights into their gene lists obtained from high-throughput experiments. one simple approach is evaluating the number of genes from a differentially expressed gene list found in a particular pathway. advanced methods calculate pathway-level statistics. in some studies, a pathway is also considered to be a gene set, a group of genes sharing common biological functions. for example, gene set enrichment analysis calculates gene set-level statistics. using gene set-level statistics has additional advantages because significant analysis at the single gene level suffers from a limited number of samples and noise. gene set based methods have also been developed to investigate phenotypic changes at the pathway level. by using pathway topology, researchers can obtain a better ranking of pathways/gene sets from their gene lists. several approaches have been proposed for constructing networks of pathways/gene sets to study complex molecular networks. yong, et al. constructed a global pathway crosstalk network and linkage network for yeast. dikla, et al. proposed a method for gleaning patterns of interactions among biological processes by analyzing protein-protein interactions, transcriptional co-expressions, and genetic interactions. gene set networks can be also applied to understanding diseases. liu, et al. proposed an approach to detect the crosstalk among alzheimer's disease related pathways and the dysfunctions in the six brain regions of ad patients. for databases of pathway/gene set network, sudhir, et al. developed human pathway database to enable the study of human pathway networks. huang, et al. developed pathway and gene-set enrichment database to further study human gene set networks. another study proposed an approach to examine interactions between pathways in mice by integrating different types of data. recently, jignesh, et al. proposed methods to construct multi-edge gene set networks to reveal insights into global relationships among biological themes or gene sets. a multi-edge gsn consists of three types of edges: co-membership edge, linkage edge, and co-enrichment edge. to our knowledge, none of the existing gene set or pathway networks contain high-level directionality information, which is essential for building models to understand the mechanistic linkages among many gene sets or pathways that are enriched in a biological condition. in this paper, we present a method of constructing a directed, regulatory gene set network called regulatory gene set network. a directed edge in an r-gsn represents a regulatory relationship from one gene set to another. compared to a linkage gene set network, which is constructed from protein interactions discarding directionality of the interactions, an r-gsn provides higher resolution of knowledge. our hypothesis is that r-gsns can reveal novel gene set relationships and provide complementary knowledge to the existing types of gsns, such as m-gsns, l-gsns, and e-gsns. in an r-gsn, a pair of gene sets are connected if a significant number of gene regulations exist between the unique genes of the gene sets. gene set and gene regulation data were collected from multiple sources. to evaluate r-gsns, m-gsns are constructed separately as baselines for each collection including kegg, reactome, and three types of gene ontology terms. we exploited the directionality information provided by r-gsns to search for significant gene sets in the network. we also compared our r-gsns with e-gsns obtained from jignesh et al. we chose alzheimer's disease as an example to study disease-specific r-gsns. in order to construct an ad specific r-gsn, gene sets from different collections were first combined into one single collection and then ad related gene sets were used to construct the ad specific r-gsn. methods in a gene set network, a node represents a gene set and an edge represents a relationship between two gene sets. a gsn helps explain biological complexity by revealing high level relationships among biological processes. we constructed a new type of gsn, a regulatory gene set network, using public available gene regulation data. another type of gsn we constructed is the co-membership gene set network, which can be a baseline network as it is constructed from annotated gene sets and thus provides experimental validation. finally, we used hypergeometric distribution to calculate significance values for each edge in both r-gsns and m-gsns. data sources the two major types of data used to construct r-gsns are gene set data and gene regulation data. gene set data provides information for nodes in a gene set network as well as information for constructing co-membership gene set network. five collections of gene sets were collected. they included kegg, reactome, go biological process, go cell component, and go molecular function from kegg and msigdb. the total number of gene sets from the five gene set collections is and the total number of genes is. gene regulation data was used for constructing r-gsns. human gene regulations from string, transfac, tred, and spike were collected. different types of id of genes obtained from different data sources were mapped to ncbi official gene symbols. in order to select only the high quality human gene regulation data, different criteria were used to filter gene regulations in different data sources. for string, gene regulations with score >= were collected; for transfac, binding site quality <= was the criterion; and for tred, only gene regulations not obtained from computational predicted method were collected. all gene regulations provided by spike were collected since they are from pathway. the total count of unique gene regulations after combining data from the four data sources is. construction of gene set networks an m-gsn was constructed separately for each gene set collection. the hypergeometric distribution was used to calculate the significance value for each co-membership edge. count the number of genes inside each of the two gene sets, gs and gs and the number of shared genes between gs and gs calculate the p-value by using the hypergeometric distribution: p-value=kkn-kn-k/nn where n is the total number of genes; n is the number of genes in gs1; k is the number of genes in gs2; and k is the number of shared genes. adjust the p-value for multiple hypotheses to control the false discovery rate by using the benjamini-hochberg procedure with p-value ≤. connect a pair of gene sets with an edge if the edge is rejected by the benjamini-hochberg procedure and thus considered as a significant edge. an r-gsn was constructed as follows: count the number of gene regulations where genes in gs regulate genes outside gs and the number of gene regulations where genes outside gs regulate genes in gs remove shared genes from gs and gs and count the number of gene regulations where the remaining genes in gs regulate the remaining genes in gs calculate the p-value by using the hypergeometric distribution using formula, where n is the total number of gene regulations; n is the number of gene regulations where genes in gs regulate genes outside gs2; k is the number of gene regulations where genes in gs are regulated from genes outside gs2; and k is the number of gene regulations from genes in gs to gs adjust the p-value for multiple hypotheses to control the false discovery rate by using the benjamini-hochberg procedure with p-value ≤. connect a pair of gene sets with a directed edge pointing from gs to gs if the edge is rejected by the benjamini-hochberg procedure and thus is significant. r-gsns were constructed separately for each gene set collection. note that it is possible that r-gsns contain a loop when both incoming and outgoing edges of a pair of gene sets are significant. r-gsns are directed networks, while m-gsns are undirected networks. before comparing r-gsns with m-gsn, the r-gsn was converted to an undirected network by discarding the directions of edges and removing loops. a disease specific gene set network a disease specific gene set network was constructed in order to show how r-gsns can help researchers explain disease complexity. alzheimer's disease was chosen to be a case study. an ad gene list containing genes from alzgene database was first collected. the five gene set collections were combined into a single global collection without removing and changing the original gene sets. from the list, ad related gene sets were selected from the global gene set collection. then the ad specific gene set network was constructed. in the network, each node represents an ad related gene set. in order to find ad related gene set, the number of genes in each gene set found in the ad gene list was counted. each ad gene list was treated as a new gene set. the same method used in constructing m-gsns was used to calculate the p-value for each gene set. there were gene sets which shared significantly high numbers of genes with the ad gene list. these ad related gene sets were collected and used to construct the ad specific r-gsns. network analysis three types of centrality were calculated for each gene set in both r-gsns and m-gsns, degree centrality, betweeness centrality, and closeness centrality. igraph software package for r was used in order to compute all network values. degree centrality represents the number of edges upon a node. a gene set with a high degree centrality is likely to be an important gene set because it acts like a hub in the network. it can also be used for comparing two different types of networks. to conduct such comparison, we calculated pearson's correlation coefficient of nodes' degree centrality between an m-gsn and an r-gsn. betweenness centrality is defined as the number of times a node acts as a bridge along the shortest path between two other nodes. in a gene set network, a gene set with a high betweenness value is likely to be a part of several biological critical paths. closeness centrality is defined as the inverse of the average length of the shortest paths between a node and all other nodes in a network. a node with a high closeness value is more central. in a gene set network, if a gene set with a high closeness value is disturbed, it is likely that a high number of gene sets will be affected by the gene set.  interpreting co-membership gene set network and regulatory gene set network when a pair of gene sets in an m-gsn are connected, the interpretation depends on gene set data types. for pathway gene sets, an m-gsn represents pathway crosstalk; whereas for go gene sets, an m-gsn represents protein moonlighting or gene sharing. when a pair of gene sets in an r-gsn are connected, they are connected by a directed edge. a directed edge in an r-gsn presents a possibility of one gene set regulating another. according to table the numbers of unique genes from all gene set collections are not much different, where reactome and go biological process contain more genes than the others. the total count of unique genes from all collections is. we calculated the proportion of regulatory edges to nodes and normalized the proportion value by using number of unique genes. we used rigi/g/si, where ri is the number of gene is sets in collection i; gi is the number of unique genes in collection i; g is the total number of unique genes; and si is the number of gene sets in collection i, to calculate the normalized proportion. in the r-gsn, go biological process has the highest normalized proportion of edges to nodes among the five gene set collections. these results indicated that pairs of biological process gene sets are more likely to have a regulatory relationship. in addition, according to the distribution of gene set sizes, the five gene set collections have similar distributions. the distributions show that the five collections contain more small gene sets than large gene sets. afor regulatory edges, a pair of gene set can have a loop if both incoming and outgoing edges are significant. the number of regulatory edges is always greater than or equal to the number of regulatory relationships because one pair of gene sets either have or do not have regulatory relationship. comparing the kegg regulatory gene set network and the kegg co-membership gene set network in the m-gsn of kegg, an edge between pathways can be considered as a pathway crosstalk. therefore, m-gsns can be baselines to compare with r-gsns. a regulatory relationship between pathways can provide complementary knowledge such as dysfunction of one pathway affecting function of other pathways. to investigate this knowledge, r-gsns and m-gsns of the kegg pathways were compared. the percentage of shared edges between m-gsns and r-gsns is the highest for kegg pathway = %). the low percentage of shared edges indicates that the r-gsn provides complementary information to the m-gsn. in addition, it is important to note that r-gsns are constructed from gene regulations which were collected from high coverage data sources. therefore, it is unlikely that r-gsns depend on the quantity and the quality of experimental data. considering both the r-gsn and m-gsn of kegg pathway gene sets, the most significant edge of the kegg r-gsn is a regulatory relationship from "cell cycle" to "cytokine-cytokine receptor interaction" with p-value = e-; whereas the co-membership edge between "cell cycle" to "cytokine-cytokine receptor interaction" has relatively low significance value. in addition, only of the top most significant regulatory edges are found in m-gsns. these findings suggest that the r-gsn reveals additional knowledge to the m-gsn. the regulatory relationship is gene set regulates gene set for the kegg r-gsn, of the most significant regulatory edges are from the "cell cycle" gene set to other kegg pathway gene sets. this suggests that changing of "cell cycle" pathway is likely to affect other pathways. this finding is corresponding to the fact that a cell cycle is a complex series of phenomena by which cellular material is duplicated and divided. therefore if a cell cycle pathway does not function appropriately, other pathways such as pathways in cancer can be affected. for the m-gsn, the kegg pathway gene sets of alzheimer's, parkinson's, and huntington's diseases have significant co-membership edges linking them together. the three co-membership edges connecting the three neurodegenerative diseases are among the top most significant co-membership edges suggesting that the three neurodegenerative diseases are highly related. in addition, of the top co-membership edges connect cancer related pathway gene sets. these edges connect the "pathways in cancer" gene set with gene sets of cancers including "small cell lung cancer", "pancreatic cancer", "melanoma", "colorectal cancer", and prostate cancer". in addition, degree centrality of each node in both networks was calculated. in the kegg m-gsn, the "pathways in cancer" gene set has the highest value of dc. in the kegg r-gsn, the "cell cycle" gene set has the highest dc and the highest out-degree centrality. the gene set which has the highest in-degree centrality is "pathways in cancer". using directionality information to calculate in-degree and out-degree centrality of the kegg r-gsn, we can find the sink, "pathways in cancer", and the source, "cell cycle", gene sets. after the dc value for each gene set was calculated, the correlation between dc of the kegg m-gsn and dc of the kegg r-gsn was calculated. the correlation coefficient is and r-squared value is. this result suggests that the gene set which is important in the m-gsn is likely to be important in the r-gsn and vice versa. we found three interesting outliers. pathway number which has regulatory dc = and co-membership dc =, is "cell cycle" suggesting that the cell cycle pathway does not tend to share genes, but to have regulatory relationship with other pathways. pathway number which has regulatory dc = and co-membership dc =, is "pathways in cancer" suggesting that pathways in cancer shared high number of genes with several pathways and their genes also regulate the unique genes of other pathways. pathway number which has regulatory dc = and co-membership dc =, is "maturity onset diabetes of the young". this pathway shares genes with only "type ii diabetes mellitus pathway". while the correlation between dc of the kegg m-gsn and dc of the kegg r-gsn is as high as, the topologies of the networks are different. this suggests that the two networks can be used to explain different phenomenon in biological systems. constructing an exclusive r-gsn we constructed a kegg exclusive r-gsn that contains only exclusive edges in the r-gsn. the correlation of dc between the kegg exclusive r-gsn and the kegg r-gsn is. the correlation of dc between the kegg exclusive r-gsn and the kegg m-gsn is, while the correlation of dc between the kegg non-exclusive r-gsn and the kegg m-gsn is. these suggest that the exclusive r-gsn can reveal important gene sets that are not likely to be revealed by m-gsn. the directionality information from a regulatory gene set network revealed "sink" and "source" gene sets in addition to "hub" gene sets, which can be revealed by constructing. table shows the top highest out-degree centrality gene sets and top highest in-degree centrality gene sets of kegg exclusive r-gsn. "cell cycle", "p signaling pathway", and "tgf-beta signalling pathway" are among the top highest out-degree centrality which can be the sources. the top highest in-degree centrality gene sets, "cell cycle", "hematopoietic cell lineage", and "cytokine-cytokine receptor interaction", are the sinks. comparison of the kegg co-enrichment network and the kegg regulatory network in the previous analysis, we used the kegg m-gsn as a baseline and compared it with the kegg r-gsn. several relationships between gene sets were found exclusively in the r-gsn. in order to validate both the kegg m-gsn and the kegg r-gsn, the two gsns were compared with the kegg co-enrichment network obtained from the study of jignesh, et al.. to construct an e-gsn, they integrate experimental gene lists and link two gene sets if the unique genes of the two gene sets are consistently enriched together across many experimentally derived gene lists. therefore, edges found in the kegg r-gsn should also be found in the kegg e-gsn. the total number of edges found in both the e-gsn and the r-gsn is, which accounts for % of the total number of edges in the e-gsn. the total number of edges found in both the e-gsn and m-gsn is which is equal to % of the total number of edges in the e-gsn. in order to calculate the significance value of the number of shared edges, the kegg e-gsn was compared with random networks. we randomly generated networks using all the gene sets from kegg. to calculate the significance value of the number of shared edges between the kegg e-gsn and the kegg m-gsn, each of the random networks contains edges, which are equal to the number of edges found in the kegg m-gsn. then fisher's exact test was used for calculating the p-value for the number of shared edges. the p-value is < e-. to calculate the significance value of the number of shared edges between the kegg e-gsn and the converted kegg r-gsn, each of the random networks contains edges. then fisher's exact test was used for calculating the p-value for the number of shared edges. the p-value is < e-. the number of shared edges between the kegg e-gsn and the kegg r-gsn is significantly high. this is corresponding to the fact that a pair of gene sets with strong regulatory relationship should be connected with a co-enrichment edge. the number of shared edges between the kegg e-gsn and the kegg m-gsn is also significantly high. this is also corresponding to the fact that a pair of gene sets with a high number of shared genes should be connected with a co-enrichment edge. a disease specific regulatory gene set network the number of shared genes between the ad gene list and each gene set in the global gene set collection was counted. out of the gene sets, have significant number of shared genes. among the ad gene sets, are from kegg; are from reactome; are from go molecular function; are from go biological process; and are from go cellular component. figure shows the degree distribution of ad-specific r-gsn constructed based on these ad gene sets. in addition, table shows the top ad related gene sets arranged by the significant values. "genes involved in lipid digestion, mobilization, and transport" gene set has the highest significant value, e- the significant value was calculated by using formula where there are genes and genes were found in the ad gene list. the ad related r-gsn of the top gene sets are presented in figure we then investigated the top degree centrality gene sets. table shows that "signal transduction" from go biological process has the highest value of dc suggesting that the signal transduction process is very important in ad. searching on pubmed found more than publications discussing the relationship between signal transduction abnormality and alzheimer's disease. considering the directionality of the ad specific r-gsn, "signal transduction" has the highest in-dc. gene sets which have the highest out-dc are "cellular protein metabolic process" and "cellular macromolecule metabolic process". these results show that directionality information from the r-gsn enables us to identify, in ad context, a sink gene set, "signal transduction", and source gene sets, "cellular protein metabolic process" and "cellular macromolecule metabolic process". furthermore, the closeness and betweenness of each gene set in the ad specific r-gsn were computed. note that the directionality of each edge in the network was considered when we calculated the shortest path for the network in order to calculate the closeness and the betweenness. for the closeness centrality, "signal transduction" from the go biological process, "pathways in cancer" from kegg, and "protein metabolic process" from the go biological process have the highest closeness values. these results suggest that inappropriate functions of the three gene sets are likely to affect high number of gene sets in ad. for the betweenness centrality, "pathways in cancer" from kegg, "system development" from go biological process, and "leishmania infection" from kegg have the highest betweenness values. these results suggest that the three gene sets are likely to be on the critical path of biological functioning for ad patients.  co-membership gene set networks and regulatory gene set networks for the five different gene set collections were constructed, compared, and studied for their biological relevance in this study. the results show that new r-gsns can provide complementary biological information to conventional m-gsns. the results also show that while the correlation between the degree centrality of the kegg m-gsn and the degree centrality of the kegg r-gsn is relatively high, the topologies of the networks are totally different. this suggests that the two networks can be employed to explain different phenomenon in biological systems. validating kegg m-gsns and kegg r-gsns separately against the kegg co-enrichment networks reveals that the numbers of shared edges are significantly high. the results suggest that the use of e-gsn to validate m-gsn and r-gsn could be a better way to describe experiments. moreover, the r-gsn specific to alzheimer's disease shows that the network can be used as a good high-level mechanistic model complementary to a more difficult-to-derive gene regulatory network model towards a systematic understanding of the disease mechanism. in this study, we describe methods for constructing both global and disease specific r-gsns. they enable future research on systems biology and translational bioinformatics. since the underlying gene regulation data are collected from high quality and high coverage data sources, the directed edges in the r-gsn do not tend to depend on the number and the quality of experimental data. moreover, directionality information from the r-gsn enables the finding of source gene sets and sink gene sets, which might be important for drug discovery or drug repositioning. tissue-specific and disease-specific gene regulations were not used in this study. therefore, the gsns obtained might be generally applicable but should be carefully reviewed and curated when they are to be used in specific disease contexts where differential gene regulations exist extensively. however, this limitation can be alleviated by constructing only condition-specific gsns instead of global gsns in the future. the framework in this study can also be later extended by collecting higher resolution such as tissue specific and disease specific gene regulation data and gene set data. separately, we have been developing a comprehensive database populated with new gene sets and all relationships identified or integrated between genes and gene sets into an online resource called pager to put the method developed in this work into a web application for biological users to explore. list of abbreviations used gsn: gene set network; m-gsn: co-membership gene set network; e-gsn: co-enrichment gene set network; l-gsn: linkage gene set network; r-gsn: regulatory gene set network; dc: degree centrality; m: co-membership edge; e: co-enrichment edge; r: regulatory edge; l: linkage edge. competing interests the authors declare that they have no competing interests. authors' contributions jake y chen conceived of this work, guided the research team by providing ideas and feedback along the way, and revised the manuscript. chayaporn collected the data, constructed gene set networks, analyzed the networks, and wrote the manuscript. liugen assisted network visualization and revised the manuscript. supplementary material additional file 1 distribution of the sizes of gene sets in each collection. this file contains five histograms presenting distribution of the sizes of gene sets in each of the five gene set collections. click here for file additional file 2 ad gene list. this file contains a list of ad related gens. click here for file additional file 3 centrality values of each gene set in the ad specific gene set network. this file contains centrality values, including degree centrality, closeness centrality, and betweenness centrality, of ad related gene sets in the ad specific r-gsn. click here for file acknowledgements we thank dr. xiaogang wu for his generous help in conceptualizing the work, providing occasional comments, and proofreading the manuscripts. we thank the support from indiana university information technology services for their partial support of this work through high-performance computing hardware and software. we thank indiana center for systems biology and personalized medicine for providing research environments that led to the timely completion of this work. declarations publication costs for this article were funded by the corresponding author's affiliated institutions, indiana university school of informatics and computing in the united states and wenzhou medical university in china. this article has been published as part of bmc genomics volume supplement 2015: selected articles from the fourth ieee international conference on computational advances in bio and medical sciences: genomics. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/16/s 
when integrated properly to functional genomics data, r-gsns can help enable future research on systems biology and translational bioinformatics. gene set networksgene regulationregulatory gene set networkspathway crosstalk2- june fourth ieee international conference on computational advances in bio and medical sciences miami beach, fl, usa several types of gene set networks have been used to study complex molecular networks such as co-membership gene set networks and co-enrichment gene set networks. we collected several gene set collections and high-quality gene regulation data in order to construct r-gsns in a comparative study with co-membership gene set networks. we described a method for constructing both global and disease-specific r-gsns and determining their significance. in this study, we proposed a new approach for constructing directed, regulatory gene set networks to reveal novel relationships among gene sets or pathways.
 traditional approaches to studying molecular networks are based on linking genes or proteins. higher-level networks linking gene sets or pathways have been proposed recently. several types of gene set networks have been used to study complex molecular networks such as co-membership gene set networks and co-enrichment gene set networks. gene set networks are useful for studying biological mechanism of diseases and drug perturbations.  in this study, we proposed a new approach for constructing directed, regulatory gene set networks to reveal novel relationships among gene sets or pathways. we collected several gene set collections and high-quality gene regulation data in order to construct r-gsns in a comparative study with co-membership gene set networks. we described a method for constructing both global and disease-specific r-gsns and determining their significance. to demonstrate the potential applications to disease biology studies, we constructed and analysed an r-gsn specifically built for alzheimer's disease.  r-gsns can provide new biological insights complementary to those derived at the protein regulatory network level or m-gsns. when integrated properly to functional genomics data, r-gsns can help enable future research on systems biology and translational bioinformatics. gene set networksgene regulationregulatory gene set networkspathway crosstalk2- june fourth ieee international conference on computational advances in bio and medical sciences miami beach, fl, usa 