BACKGROUND
from the perspective of network medicine, a disease phenotype is rarely a consequence of an abnormality in a single biomolecule , but reflects various pathobiological processes that interact in a complex network  <cit> .
one single factor can exert certain effects on disease when studying it alone, while this effect may be vanished when studying it within one network or pathway  <cit> , and vice versa.
therefore, biomolecules should be studied in the context of biological systems they are involved in  <cit> .
perhaps the abstraction for a biological system is network, such as transcriptional regulatory networks, signal transduction networks, protein interaction networks and metabolic networks  <cit> .
in the biological networks, the vertices represent biomolecules, and the edges represent functional, causal or physical interactions between the vertices.
different types of networks are often used to represent diverse types of biological processes, each of which stores information about levels and interactions related to specific biomolecules  <cit> .
in fact, different physiological conditions may manifest as different networks.
moreover, complex disease are multi-factorial and analyzing the individual components is insufficient, so it is essential to dissect how these components interact with each other and weave into one network, and how these interactions differ with respect to disease status.
statistical comparison of group difference in biological networks or pathways can provide new insight into the underlying disease mechanism, and have extensive biomedical and clinical applications .
for instance, a better understanding of the effects of molecular interconnectedness on disease progression may lead to superior identification of disease related biomolecules and pathways, which may further offer more effective targets for drug development in a cost-effective and timely manner.
on the other hand, identifying biological and environmental causes of human diseases has always been one of the central concerns in epidemiology.
however, traditional epidemiology has been pejoratively labeled as the “black box” epidemiology  <cit> , and increasingly suffered from criticism partly due to the fact that too much attention has been paid to the identification of a single risk factor rather than the network or pathway related to a disease, which led to difficulty to deeply explore disease mechanism  <cit> .
it is highly desirable to unlock the black box underlying observed associations and to illuminate the biological interaction mechanisms of disease-related components hiding behind the black box.
there are unmet needs to access multi-level omics data on the population level.
thanks to the development of recent technological advances in high-throughput omics platforms, we can enable the acquisition of omics data at unprecedented speed and amounts, and further integrate various omics data with traditional epidemiology to promote the development of systems epidemiology  <cit> .
it offers the potential to provide new insight into the underlying disease mechanisms in breadth and depth at human population level.
under the framework of systems epidemiology, the focus has been shifted from identification of single factor to exploration of specific networks or pathways contributing to disease  <cit> .
in a word, it is in great needs to do statistical comparison of biological networks.
so far, several methods have been proposed to utilize network topology information to carry out various biomedical tasks.
langfelder et al.
<cit>  provided several measures for comparing network topologies for weighted correlation networks.
zhang et al.
<cit>  proposed a differential dependency network analysis to detect topological changes in transcriptional networks between subclasses of breast cancer.
valcarcel et al.
<cit>  introduced a formal statistical method for the differential analysis of molecular pair-wise associations via network representation.
recently, yates et al.
<cit>  developed an additive element-wise-based dissimilarity measure for biological network hypothesis tests.
however, most of above methods mainly focus on the difference of network topology and are unable to account for the changes of vertices.
although in most situations, the differences of single vertices-wise or edges-wise may be weak, their aggregated differences can be quite strong.
it will undoubtedly lose statistical power to only consider the connection with the topological difference between two networks.
meanwhile, non-parametric permutation procedures are commonly employed to perform analysis in most existed methods, which were inevitably time-consuming, especially for big data.
the premise for networks or pathways comparison is to make clear the cause of biological network difference.
generally, both changes in the nodes level , and changes in the edges  can lead to the whole network difference.
reverter et al.
<cit>  presented an analytical procedure to simultaneously identify genes that were differentially expressed  as well as genes that are differentially connected  for unweighted networks.
their methods depend heavily on the specific definition of de and dc, and the two-component mixture of bi-variate normal distribution may be violated in other biological networks, though it may be reasonable in gene expression network.
furthermore, weighted  networks are commonly encountered and increasingly relevant in biological applications .
statistical methods for detecting the group difference in weighted biological networks are still in great demand.
in this article, we proposed a new score-based network difference measure  as a powerful test statistic to detect group difference in weighted networks, which simultaneously capture the difference of vertices and edges.
various simulations were conducted to evaluate its type i error and statistical power, compared with other existed method.
two real data sets about gwas of leprosy and gene expression of ovarian cancer were further analyzed to show their performance in practice.
methods
statistical model
a weighted biological network can be modeled as an undirected graph g = , where v is the set of vertices  and e is the set of edges .
two vertices, representing biomolecules, are connected by an undirected edge if there is an association between them.
each edge can be assigned a weight resembling the strength of evidence for the association.
we denote the two networks in two groups  by gd and gc respectively, suppose both gd and gc have the same number of vertices  and edges , the null hypothesis test is h0 : gd = gc.
let v and e denote the set of all vertices and edges in gd, xidxjd indicate the edge xid−xjd , βijd represent the strength of association between xid and xjd if xidxjd existed.
for individual l , the trait value is denoted as yl, yl= <dig> l∈case <dig> l∈control and the ith vertex is denoted as xli.
under h <dig>  networks in two groups are identical not only in the average vertices levels but also in the connection strength.
the score test vector of vertices is dv = t, where div=∑l=1nyl−y¯xli measures the contribution of vertex xi to the disease.
analogously, the score test vector of edges is de = t, where dke=∑l=1nyl−y¯xli−x¯ixlj−x¯j measures the contribution of connection strength between xi and xj  to the disease.
then the proposed overall network difference measure can be defined as netdifm=dtΣ−1d, where d=dvde,Σ = cov =  × , p, q =  <dig>   <dig>  ⋯, .
the estimated covariance matrix of d can be represented as  and calculated as follows,
1) for Σv, p, q =  <dig>   <dig>  ⋯, m, σpq=∑l=1nyl−y¯2covxpxq,xp=x1px2p⋯xnp;
2) for Σe, p, q = m +  <dig>  m +  <dig>  ⋯, m + k, σpq=∑l=1nyl−y¯2covzpzq,zp=xi−x¯i×xj−x¯j;
3) for Σve , p =  <dig>   <dig>  ⋯, m , q = m +  <dig>  m +  <dig>  ⋯, m + k σpq=∑l=1nyl−y¯2covxpzq.
naturally, for a large sample size, netdifm has a centered χ <dig> distribution under the null hypothesis .
when sample size is small, a permutation procedure can be conducted as follows to get the empirical p value and assess the statistical significance.
calculate the test statistic netdifm from the original sample;  randomly assign subjects to one of two groups, the sample size for each group keeps the same as the original data;  perform the above steps q times and calculate the test statistic for each repeated sample, netdifmi*, i =  <dig>   <dig>  ⋯, q;  estimate the p value according to p−value=1q∑i=1qinetdifmi*>netdifm, where i is the indicator function.
intuitively, considering the elements of one network are not more than vertices and edges, an element-wise measure may be expected to have the ability to identify the group difference in biological networks.
a vertices and edges wise difference measure , through the simple summation of vertices difference and edges, can be constructed as vewdm=1m∑i=1mti2+1k∑i=1m∑j≠imuij <dig> where ti=x¯id−x¯icvarx¯id+varx¯ic, x¯id and x¯ic indicate the sample mean of xi in gd and gc respectively; uij=βijd−βijcvarβijd+varβijc, when the strength of edges are quantified by the pearson correlations rij, uij=zijd−zijc/1nd−3+1nc− <dig>  zij are the fisher-transforms of the correlations zij=12ln1+rij1−rij, nd and nc are the corresponding sample sizes.
the proposed vewdm seems to be the linear combination of some chi-square statistics.
actually, the asymptotic theoretical properties have been explored for the linear combination of independent chi-squares  <cit> .
nevertheless, it is quite complex and difficult here to obtain the asymptotic distribution of vewdm, since the correlations between different vertices and different edges statistics  highly depend on the specific network structure.
in other words, the asymptotic properties might be network-specific.
to solve this problem, we adopted the strategy of a permutation test to make statistical inference.
simulation
simulations were designed to evaluate the type i error rate and statistical power, and to compare the performance of netdifm, vewdm and yates’d  under different sample size and network topological structure.
the statistical power is defined as the probability that the test correctly rejects the null hypothesis  when the alternative hypothesis  is true.
it can be estimated from the empirical distribution as the proportion of observations for which the p-value is less than given nominal level .
for the specific network with m vertices and k edges, the simulated m-dimensional variables  were generated from a multivariate normal distribution nm with mean vector μ and covariance matrix Σ using the r “mvtnorm” package.
we specified the mean vector μ =  and covariance matrix Σ = m × m, where iij= <dig> xixj∈eg <dig> xixj∉eg, i ≠ j, i, j =  <dig>   <dig>  ⋯, m was the indicator function.
under the null hypotheses , the data was generated by setting μd = μc and iijdβijd = iijcβijc.
<dig> simulations were repeated to assess the type i error of the above methods given various sample sizes under different network scale, including network with ten vertices and  <dig> edges  and another one with  <dig> vertices and  <dig> edges .
under the alternative hypotheses, three scenarios were considered.fig.
<dig> weighted biological networks.
a an assumed network including ten vertices and  <dig> edges.
b an assumed network including  <dig> vertices and  <dig> edges
for scenario  <dig>  only vertices  were different between gd and gc.
we simulated two vertices difference with μ3d − μ3c =  <dig>  and μ10d − μ10c =  <dig>  under the network topological structure as in fig.
1a.
three vertices difference with μ6d − μ6c =  <dig> , μ8d − μ8c =  <dig>  and μ17d − μ17c =  <dig>  were also designed under the relative larger scale network as in fig.
1b.
for scenario  <dig>  only edges  were different between gd and gc.
we simulated three edges difference with β35d − β35c = −  <dig> , β57d − β57c =  <dig>  and β <dig> d − β <dig> c =  <dig>  under the network topological structure as in fig.
1a.
seven edges difference with β <dig> c − β <dig> d = −  <dig> , β <dig> c − β <dig> d = β <dig> c − β <dig> d = β <dig> c − β <dig> d = β <dig> c − β <dig> d =  <dig>  and β23c − β23d = β <dig> c − β <dig> d =  <dig>  were also designed under the relative larger scale network as in fig.
1b.
for scenario  <dig>  both vertices and edges were designed to be different between gd and gc.
under the topology structure as in fig.
1a, we combined the settings in scenario  <dig> and scenario  <dig> , so as for the topology structure as in fig.
1b.
for each scenario,  <dig> replicates were used to evaluate statistical power.
p-values of the proposed netdifm were assessed using both the asymptotic distribution and the empirical null distribution obtained from  <dig> times permutations.
it is necessary to assess the performance of the proposed statistics, given the deviation from the normal distribution.
for the network with ten vertices and  <dig> edges, we designed the following two scenarios,  conduct the exponential transformation for five vertices randomly chosen among the ten vertices;  do the exponential transformation for all ten vertices.
for each scenario, we evaluate the type i error rate and statistical power under the same three scenarios mentioned as above.
furthermore, we also provided estimated computing time under different network with sample size  <dig> and  <dig> permutations, using one laptop as an intel pentiumt <dig> with a  <dig>  ghz cpu and  <dig> gb ram.
application
gwas data of leprosy
by ingenuity pathways analysis, a plausible biologic network underlying susceptibility to leprosy was provided for depicting the functional relationship between some susceptibility genes identified from gwas of leprosy  <cit> .
using the initial gwas data with  <dig> cases and  <dig> controls, we attempted to detect the difference of the networks including genes card <dig>  hla-drb <dig>  ripk <dig>  card <dig> and ifng.
all participants provided written informed consent, and the study was approved by the ethics committees of shandong academy of medical science  <cit> .
these five genes located on different chromosomes and totally contained  <dig> snps , with network structure given in fig.
2a.
since each gene contained several snps, we first employed principal component analysis and conducted the statistical network comparison by treating the first principal component as the network vertices.fig.
<dig> gene networks.
a gene interaction network of leprosy.
b subnetwork of pi3k-aktsignaling pathway.
c subnetwork of notch signaling pathway
gene expression data of ovarian cancer
tothillet al.
<cit>  used high-density expression oligonucleotide microarrays for profiling  <dig> well-annotated serous and endometrioid invasive ovarian, fallopian tube, and peritoneal cancers.
the subjects were divided into a c <dig> subtype, with  <dig> patients, and a c2–c <dig> subtype, with  <dig> patients.
complete expression data are available on geo .
the proposed method was also applied to detect the network difference between these two groups .
here we studied two specific pathways  reported in the literatures  to be relevant to ovarian cancer.
the subnetwork of pi3k-akt signaling pathway from the kegg pathway database with  <dig> genes contained  <dig> probe sets  was abstracted into network with topological structure shown in fig.
2b.
the subnetwork of notch signaling pathway with  <dig> genes contained  <dig> probe sets  was abstracted into network with topological structure shown in fig.
2c.
all probe sets corresponding to the same gene symbol were first averaged to obtain gene-level expression measurements.
RESULTS
simulation
shown in table  <dig> are the estimated type i error rates of the proposed netdifm, vewdm, netdifm based on permutation  and yates’d under different sample sizes.
it reveals that all type i error rates based on permutation procedure are close to given nominal level .
netdifm tended to be slightly conservative under small sample size, while using the asymptotic distribution maintains a good control of type i error rate under large sample size.table  <dig> type i error rates of four methods
netdifm
vewdm
yates’d
figure 3a indicates the statistical power under scenario  <dig> when only vertices changed with the network topological structure demonstrates in fig.
1a.
as expected, yates’d has no power due to that it can only capture the edge change.
netdifm is substantially more powerful than vewdm, and it is slightly less powerful than its permutation-based type.
similar trend could also be found under the relative larger scale network .fig.
<dig> the statistical power of the four methods under the scenario only vertices was different between two groups.
a the power under the network topological structure as in fig.
1a.
b the power under the network topological structure as in fig.
1b
shown in fig.
<dig> is the performance under scenario  <dig> .
the statistical power of all methods monotonically increases with sample size.
netdifm has much higher power than that of vewdm and yates’d.
the power of netdifm and yates’d keep almost the same in the larger scale network .fig.
<dig> the statistical power of the four methods under the scenario only edges was different between two groups.
a the power under the network topological structure as in fig.
1a.
b the power under the network topological structure as fig.
1b
figure  <dig> illustrates the statistical power under the scenario  <dig> .
both netdifm and vewdm are much more powerful than yates’d, and netdifm still has the best performance.fig.
<dig> the statistical power of the four methods under the scenario both vertices and edges were different between two groups.
a the power under the network topological structure as in fig.
1a.
b the power under the network topological structure as in fig.
1b
to evaluate the scalability and computational efficiency of the proposed methods, we also conducted simulations using a larger network with  <dig> vertices and  <dig> edges .
it is clear that the proposed netdifm still have the best performance .
figure  <dig> indicates the results given the deviation from the normal distribution, where the proposed statistics still hold the relative better performance than other method.fig.
<dig> the statistical power of the four methods given the deviation from the normal distribution.
a, b, c the power under the scenario conducting the exponential transformation for five vertices, when only vertices change , only edge changes , both vertices and edges change .
d, e, f the power under the scenario doing the exponential transformation for all ten vertices, when only vertices change , only edge changes , both vertices and edges change
table  <dig> presents the estimated computing time.
it indicates that the proposed netdifm indeed runs fast, and the computational time increases as the network become larger.table  <dig> computing time  of four methods with sample size  <dig> and  <dig> permutations
netdifm
vewdm
yates’d
network <dig> ten vertices and  <dig> edges; network <dig>  <dig> vertices and  <dig> edges;
network <dig>  <dig> vertices and  <dig> edges
the results of application
network difference analysis for both the gwas of leprosy and gene expression data of ovarian cancer further confirm in practice that the proposed netdifm captured the network changes.
shown in table  <dig> are the results of the proposed netdifm and other methods for detecting the network difference between two groups.
the difference of gene interaction network with  <dig> genes can be detected significantly at α =  <dig>  by netdifm, netdifmpm and vewdm.table  <dig>
p-values of the four methods for the two real data sets
netdifm
vewdm
yates’d
group difference of the subnetwork of pi3k-akt signaling pathway was detected significantly at α =  <dig>  by netdifm, netdifmpm and vewdm.
when applied to the subnetwork of notch signaling pathway, all four methods can detect the network difference significantly .
shown in table  <dig> are the results of the proposed netdifm and other methods for detecting the specific vertices, treating a vertex as well as its connected edge as a network, under <dig> permutation times.table  <dig>
p-values of the four methods for detecting the genes in the pi3k-akt signaling pathway and notch signaling pathway
netdifm
vewdm
yates’d
kras
pik3ca
irs1
pdpk1
akt1
mdm2
tp53
dll3
dtx2
crebbp
ptcra
jag1
dvl2
snw1
hes1
rbpj
notch2
psenen
adam17
numb
ncor2
discussion and
CONCLUSIONS
complex disease is largely determined by a number of biomolecules interwoven into networks, rather than a single biomolecule.
group-level comparison of network properties  may shed light on underlying biological processes or disease mechanisms, and benefit the design of drug targets and drug combination for the therapy of complex diseases.
meanwhile, although the conventional single-based paradigm has successfully identified a list of risk factors, one common sense is that there still exist an intermediate “black box” between the exposures and the disease phenotypes .
in the “black box”, various risk factors weaved into complicated biological networks dominating the disease occurrence, development and prognosis.
recent advances in high-throughput technologies and omics resources are revolutionizing biomedical research, and allow a transition from the traditional paradigm for biological and epidemiological studies of complex diseases to a new paradigm based on systems epidemiology  <cit> .
under this framework, network-based analysis has been integrated into observational study designs to organize the interdependencies of biomolecules and other factors at a human population level, expecting to open the “black box”.
a key but inadequately addressed issue is still to develop valid statistical method to test possible differences of the networks between two groups.
in our previous study  <cit> , we have developed a statistical method for detecting the pathway effect contributing to disease, mainly under the framework of systems epidemiology.
however, this method is limited to the pathway with chain structure, and can only capture the edge changes while omitting the vertex changes.
meanwhile, the nonparametric bootstrap method has to be used to obtain the significance.
at present study, we proposed a score-based powerful statistical test to detect the significant changes in biological networks between two different conditions .
it can simultaneously capture the vertex changes and edge changes.
various simulations were conducted to assess the reliability and statistical power of the proposed method.
it indicated that both netdifm and vewdm were much more powerful than yates’d, and netdifm kept the best performance under various scenarios , and it can indeed capture the perturbation of vertices and edges in the network simultaneously.
one strength for netdifm is that we can obtain its theoretic property, and thus can avoid the high computation burden.
as expected, the proposed netdifm indeed runs fast .
the vewdm was used fisher r-to-z transformation to identify significant differences between two correlations.
fukushima et al.
<cit>  also developed an r package to identify differential correlations between two conditions based on fisher’s z-test which affords users a simple and effective framework in omics data.
the vewdm can be treated as a global measure to detect the group difference of networks between two conditions, accounting for not only edges difference but also vertices difference.
even though one is interested in testing particular vertex or edge rather than the whole network, its connected edge should also be considered.
two real data sets analyses further highlighted that netdifm had more advantage in practice.
in the gwas data of leprosy, we detected a candidate gene interaction network containing five genes.
for the gene expression data of ovarian cancer, two candidate subnetworks, pi3k-akt signaling pathway and notch signaling pathway, respectively were considered and identified, suggesting that the proposed method is capable of identifying differential gene expression and gene-gene co-expression patterns, which are certainly helpful for us to further understand the underlying disease mechanism.
rao et al.
<cit>  reported that combined overexpression of ova <dig> and mdm <dig> promotes oncogenesis by enhancing activation of the igf-1r–erk1/ <dig> signaling pathway, and jag <dig> enhances ovarian cancer cell growth and cisplatin-resistance  <cit> .
the expression of hes <dig> is confirmed to be strongly associated with the pathogenesis of ovarian endometriomas  <cit> .
meanwhile, decreased notch <dig> expression is associated with the poorly differentiated serous epithelial ovarian carcinoma histology  <cit> .
rbpj underexpression in ovarian tumor tissue relative to matched normal tissue  <cit> .
moreover, adam <dig> is one of the several metalloproteinases that play a key role in epidermal growth factor receptor signalling and can be a potential target antigen to devise novel immunotherapeutic strategies against ovarian cancer  <cit> .
the pi3k-akt and notch pathways and their abundant associated genes comprise complicated networks, which play a significant role in the progressive growth of tumor cells.
network difference can result from not only changes of vertices but also changes of edges, and the changes of vertices-wise and edges-wise are often closely related.
for instance, differential expression of genes may be due to either mutation of its own gene or the effects of expression changes of other genes in the network.
however, the degree of differential expression of one gene due to its own mutation is often lower than affected by expressions of upstream genes in the network  <cit> .
reverter et al.
<cit>  presented an analytical procedure to simultaneously identify differential gene expression and connectivity for unweighted gene network.
in their work, an edge between two genes is established if the absolute value of the correlation coefficient exceeds a fixed threshold.
consequently, if we set the threshold less than  <dig> , and the correlation coefficient between gene a and gene b is  <dig>  in cases and  <dig>  in controls, then the connection between gene a and gene b is treated as no difference between cases and controls.
while in this situation, there exists a difference of the strength of connection between gene a and gene b among cases and controls, given our methods focus mainly on weighted biological networks.
furthermore, the covariance structure between vertex changes and edge changes has been embedded into the proposed score-based network difference measure.
in addition, one would be more interested in testing particular vertex or edge  rather than the whole network or pathway.
actually, a vertex as well as its connected edge can be treated as a subnetwork, and the proposed network difference method can easily be extended to identify the specific vertices.
even though some local interventions were often generated to prevent and cure a particular disease, it is essential to understand the global system.
the ‘think globally, act locally’ paradigm should be strongly embedded into our mind  <cit> .
one limitation in our paper is that we assume the network topological structure is fixed, and little attention has been paid on the network structure learning problem.
constructing network structure means determining every possible edge with highest degree of data matching, and often one joint probability distribution of a number of variables can reflect more than one network structure.
usually, combining experimental evidence with their experience, most biologists and clinical researchers have a growing awareness of the interplay between the biological components and can depict more or less the specific network or pathway for the corresponding biological process.
meanwhile, numerous databases  can be further borrowed to establish the network structure.
the proposed netdifm will do not work in its current version under the scenario when the covariance matrix is not invertible.
one possible solution is to first apply a shrinkage strategy to simplify the network, and then adopt the proposed statistic.
for instance, we can first remove those edges if the correlation between the two linked vertices is smaller than a predefined threshold, and then apply the proposed test to the remaining network.
statistical comparisons of group difference in biological networks are highly desirable.
the proposed network difference measure was valid and powerful to detect biological network difference.
r code to perform netdifm, netdifmpm and vewdm is provided in the additional file  <dig>
availability of supporting data
the gene expression data of ovarian cancer were downloaded from the geo datasets .
additional files
additional file 1: table s <dig>  the location and snp number for five susceptibility genes belonging to the network associated with leprosy.
table s <dig>  the location and probe sets number for genes belonging to the subnetwork from pi3k-akt signaling pathway and notch signaling pathway.
table s <dig>  type i error rates of four methods .
table s <dig>  type i error rates of four methods given the deviation from the normal distribution.
figure s <dig>  a network including  <dig> vertices and  <dig> edges.
figure s <dig>  the statistical power of the four methods under three scenarios.
only vertices change,  only edges change,  both vertices and edges change.
additional file 2:
r code.
abbreviations
dcdifferentially connected
dedifferentially expressed
gwasgenome-wide association study
keggkyoto encyclopedia of genes and genomes
netdifmnetwork difference measure
netdifmpmnetdifm based on permutation
vewdmvertices and edges wise difference measure
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jj, zy and fx conceived, designed the study and implemented the data analysis, jj, zy and xz drafted the manuscript.
all authors read and approved the final manuscript.

 from the perspective of network medicine, a disease phenotype is rarely a consequence of an abnormality in a single biomolecule, but reflects various pathobiological processes that interact in a complex network. one single factor can exert certain effects on disease when studying it alone, while this effect may be vanished when studying it within one network or pathway, and vice versa. therefore, biomolecules should be studied in the context of biological systems they are involved in. perhaps the abstraction for a biological system is network, such as transcriptional regulatory networks, signal transduction networks, protein interaction networks and metabolic networks. in the biological networks, the vertices represent biomolecules, and the edges represent functional, causal or physical interactions between the vertices. different types of networks are often used to represent diverse types of biological processes, each of which stores information about levels and interactions related to specific biomolecules. in fact, different physiological conditions may manifest as different networks. moreover, complex disease are multi-factorial and analyzing the individual components is insufficient, so it is essential to dissect how these components interact with each other and weave into one network, and how these interactions differ with respect to disease status. statistical comparison of group difference in biological networks or pathways can provide new insight into the underlying disease mechanism, and have extensive biomedical and clinical applications. for instance, a better understanding of the effects of molecular interconnectedness on disease progression may lead to superior identification of disease related biomolecules and pathways, which may further offer more effective targets for drug development in a cost-effective and timely manner. on the other hand, identifying biological and environmental causes of human diseases has always been one of the central concerns in epidemiology. however, traditional epidemiology has been pejoratively labeled as the “black box” epidemiology, and increasingly suffered from criticism partly due to the fact that too much attention has been paid to the identification of a single risk factor rather than the network or pathway related to a disease, which led to difficulty to deeply explore disease mechanism. it is highly desirable to unlock the black box underlying observed associations and to illuminate the biological interaction mechanisms of disease-related components hiding behind the black box. there are unmet needs to access multi-level omics data on the population level. thanks to the development of recent technological advances in high-throughput omics platforms, we can enable the acquisition of omics data at unprecedented speed and amounts, and further integrate various omics data with traditional epidemiology to promote the development of systems epidemiology. it offers the potential to provide new insight into the underlying disease mechanisms in breadth and depth at human population level. under the framework of systems epidemiology, the focus has been shifted from identification of single factor to exploration of specific networks or pathways contributing to disease. in a word, it is in great needs to do statistical comparison of biological networks. so far, several methods have been proposed to utilize network topology information to carry out various biomedical tasks. langfelder et al. provided several measures for comparing network topologies for weighted correlation networks. zhang et al. proposed a differential dependency network analysis to detect topological changes in transcriptional networks between subclasses of breast cancer. valcarcel et al. introduced a formal statistical method for the differential analysis of molecular pair-wise associations via network representation. recently, yates et al. developed an additive element-wise-based dissimilarity measure for biological network hypothesis tests. however, most of above methods mainly focus on the difference of network topology and are unable to account for the changes of vertices. although in most situations, the differences of single vertices-wise or edges-wise may be weak, their aggregated differences can be quite strong. it will undoubtedly lose statistical power to only consider the connection with the topological difference between two networks. meanwhile, non-parametric permutation procedures are commonly employed to perform analysis in most existed methods, which were inevitably time-consuming, especially for big data. the premise for networks or pathways comparison is to make clear the cause of biological network difference. generally, both changes in the nodes level, and changes in the edges can lead to the whole network difference. reverter et al. presented an analytical procedure to simultaneously identify genes that were differentially expressed as well as genes that are differentially connected for unweighted networks. their methods depend heavily on the specific definition of de and dc, and the two-component mixture of bi-variate normal distribution may be violated in other biological networks, though it may be reasonable in gene expression network. furthermore, weighted networks are commonly encountered and increasingly relevant in biological applications. statistical methods for detecting the group difference in weighted biological networks are still in great demand. in this article, we proposed a new score-based network difference measure as a powerful test statistic to detect group difference in weighted networks, which simultaneously capture the difference of vertices and edges. various simulations were conducted to evaluate its type i error and statistical power, compared with other existed method. two real data sets about gwas of leprosy and gene expression of ovarian cancer were further analyzed to show their performance in practice. methods statistical model a weighted biological network can be modeled as an undirected graph g =, where v is the set of vertices and e is the set of edges. two vertices, representing biomolecules, are connected by an undirected edge if there is an association between them. each edge can be assigned a weight resembling the strength of evidence for the association. we denote the two networks in two groups by gd and gc respectively, suppose both gd and gc have the same number of vertices and edges, the null hypothesis test is h0: gd = gc. let v and e denote the set of all vertices and edges in gd, xidxjd indicate the edge xid−xjd, βijd represent the strength of association between xid and xjd if xidxjd existed. for individual l, the trait value is denoted as yl, yl= l∈case l∈control and the ith vertex is denoted as xli. under h networks in two groups are identical not only in the average vertices levels but also in the connection strength. the score test vector of vertices is dv = t, where div=∑l=1nyl−y¯xli measures the contribution of vertex xi to the disease. analogously, the score test vector of edges is de = t, where dke=∑l=1nyl−y¯xli−x¯ixlj−x¯j measures the contribution of connection strength between xi and xj to the disease. then the proposed overall network difference measure can be defined as netdifm=dtΣ−1d, where d=dvde,Σ = cov =  ×, p, q =     ⋯,. the estimated covariance matrix of d can be represented as and calculated as follows, 1) for Σv, p, q =     ⋯, m, σpq=∑l=1nyl−y¯2covxpxq,xp=x1px2p⋯xnp; 2) for Σe, p, q = m +   m +   ⋯, m + k, σpq=∑l=1nyl−y¯2covzpzq,zp=xi−x¯i×xj−x¯j; 3) for Σve, p =     ⋯, m, q = m +   m +   ⋯, m + k σpq=∑l=1nyl−y¯2covxpzq. naturally, for a large sample size, netdifm has a centered χ distribution under the null hypothesis. when sample size is small, a permutation procedure can be conducted as follows to get the empirical p value and assess the statistical significance. calculate the test statistic netdifm from the original sample; randomly assign subjects to one of two groups, the sample size for each group keeps the same as the original data; perform the above steps q times and calculate the test statistic for each repeated sample, netdifmi*, i =     ⋯, q; estimate the p value according to p−value=1q∑i=1qinetdifmi*>netdifm, where i is the indicator function. intuitively, considering the elements of one network are not more than vertices and edges, an element-wise measure may be expected to have the ability to identify the group difference in biological networks. a vertices and edges wise difference measure, through the simple summation of vertices difference and edges, can be constructed as vewdm=1m∑i=1mti2+1k∑i=1m∑j≠imuij where ti=x¯id−x¯icvarx¯id+varx¯ic, x¯id and x¯ic indicate the sample mean of xi in gd and gc respectively; uij=βijd−βijcvarβijd+varβijc, when the strength of edges are quantified by the pearson correlations rij, uij=zijd−zijc/1nd−3+1nc− zij are the fisher-transforms of the correlations zij=12ln1+rij1−rij, nd and nc are the corresponding sample sizes. the proposed vewdm seems to be the linear combination of some chi-square statistics. actually, the asymptotic theoretical properties have been explored for the linear combination of independent chi-squares. nevertheless, it is quite complex and difficult here to obtain the asymptotic distribution of vewdm, since the correlations between different vertices and different edges statistics highly depend on the specific network structure. in other words, the asymptotic properties might be network-specific. to solve this problem, we adopted the strategy of a permutation test to make statistical inference. simulation simulations were designed to evaluate the type i error rate and statistical power, and to compare the performance of netdifm, vewdm and yates’d under different sample size and network topological structure. the statistical power is defined as the probability that the test correctly rejects the null hypothesis when the alternative hypothesis is true. it can be estimated from the empirical distribution as the proportion of observations for which the p-value is less than given nominal level. for the specific network with m vertices and k edges, the simulated m-dimensional variables were generated from a multivariate normal distribution nm with mean vector μ and covariance matrix Σ using the r “mvtnorm” package. we specified the mean vector μ =  and covariance matrix Σ = m × m, where iij= xixj∈eg xixj∉eg, i ≠ j, i, j =     ⋯, m was the indicator function. under the null hypotheses, the data was generated by setting μd = μc and iijdβijd = iijcβijc. simulations were repeated to assess the type i error of the above methods given various sample sizes under different network scale, including network with ten vertices and edges and another one with vertices and edges. under the alternative hypotheses, three scenarios were considered.fig. weighted biological networks. a an assumed network including ten vertices and edges. b an assumed network including vertices and edges for scenario only vertices were different between gd and gc. we simulated two vertices difference with μ3d − μ3c =  and μ10d − μ10c =  under the network topological structure as in fig. 1a. three vertices difference with μ6d − μ6c = , μ8d − μ8c =  and μ17d − μ17c =  were also designed under the relative larger scale network as in fig. 1b. for scenario only edges were different between gd and gc. we simulated three edges difference with β35d − β35c = − , β57d − β57c =  and β d − β c =  under the network topological structure as in fig. 1a. seven edges difference with β c − β d = − , β c − β d = β c − β d = β c − β d = β c − β d =  and β23c − β23d = β c − β d =  were also designed under the relative larger scale network as in fig. 1b. for scenario both vertices and edges were designed to be different between gd and gc. under the topology structure as in fig. 1a, we combined the settings in scenario and scenario, so as for the topology structure as in fig. 1b. for each scenario, replicates were used to evaluate statistical power. p-values of the proposed netdifm were assessed using both the asymptotic distribution and the empirical null distribution obtained from times permutations. it is necessary to assess the performance of the proposed statistics, given the deviation from the normal distribution. for the network with ten vertices and edges, we designed the following two scenarios, conduct the exponential transformation for five vertices randomly chosen among the ten vertices; do the exponential transformation for all ten vertices. for each scenario, we evaluate the type i error rate and statistical power under the same three scenarios mentioned as above. furthermore, we also provided estimated computing time under different network with sample size and permutations, using one laptop as an intel pentiumt with a ghz cpu and gb ram. application gwas data of leprosy by ingenuity pathways analysis, a plausible biologic network underlying susceptibility to leprosy was provided for depicting the functional relationship between some susceptibility genes identified from gwas of leprosy. using the initial gwas data with cases and controls, we attempted to detect the difference of the networks including genes card hla-drb ripk card and ifng. all participants provided written informed consent, and the study was approved by the ethics committees of shandong academy of medical science. these five genes located on different chromosomes and totally contained snps, with network structure given in fig. 2a. since each gene contained several snps, we first employed principal component analysis and conducted the statistical network comparison by treating the first principal component as the network vertices.fig. gene networks. a gene interaction network of leprosy. b subnetwork of pi3k-aktsignaling pathway. c subnetwork of notch signaling pathway gene expression data of ovarian cancer tothillet al. used high-density expression oligonucleotide microarrays for profiling well-annotated serous and endometrioid invasive ovarian, fallopian tube, and peritoneal cancers. the subjects were divided into a c subtype, with patients, and a c2–c subtype, with patients. complete expression data are available on geo. the proposed method was also applied to detect the network difference between these two groups. here we studied two specific pathways reported in the literatures to be relevant to ovarian cancer. the subnetwork of pi3k-akt signaling pathway from the kegg pathway database with genes contained probe sets was abstracted into network with topological structure shown in fig. 2b. the subnetwork of notch signaling pathway with genes contained probe sets was abstracted into network with topological structure shown in fig. 2c. all probe sets corresponding to the same gene symbol were first averaged to obtain gene-level expression measurements.  simulation shown in table  are the estimated type i error rates of the proposed netdifm, vewdm, netdifm based on permutation and yates’d under different sample sizes. it reveals that all type i error rates based on permutation procedure are close to given nominal level. netdifm tended to be slightly conservative under small sample size, while using the asymptotic distribution maintains a good control of type i error rate under large sample size.table type i error rates of four methods netdifm vewdm yates’d figure 3a indicates the statistical power under scenario when only vertices changed with the network topological structure demonstrates in fig. 1a. as expected, yates’d has no power due to that it can only capture the edge change. netdifm is substantially more powerful than vewdm, and it is slightly less powerful than its permutation-based type. similar trend could also be found under the relative larger scale network.fig. the statistical power of the four methods under the scenario only vertices was different between two groups. a the power under the network topological structure as in fig. 1a. b the power under the network topological structure as in fig. 1b shown in fig.  is the performance under scenario. the statistical power of all methods monotonically increases with sample size. netdifm has much higher power than that of vewdm and yates’d. the power of netdifm and yates’d keep almost the same in the larger scale network.fig. the statistical power of the four methods under the scenario only edges was different between two groups. a the power under the network topological structure as in fig. 1a. b the power under the network topological structure as fig. 1b figure  illustrates the statistical power under the scenario. both netdifm and vewdm are much more powerful than yates’d, and netdifm still has the best performance.fig. the statistical power of the four methods under the scenario both vertices and edges were different between two groups. a the power under the network topological structure as in fig. 1a. b the power under the network topological structure as in fig. 1b to evaluate the scalability and computational efficiency of the proposed methods, we also conducted simulations using a larger network with vertices and edges. it is clear that the proposed netdifm still have the best performance. figure  indicates the results given the deviation from the normal distribution, where the proposed statistics still hold the relative better performance than other method.fig. the statistical power of the four methods given the deviation from the normal distribution. a, b, c the power under the scenario conducting the exponential transformation for five vertices, when only vertices change, only edge changes, both vertices and edges change. d, e, f the power under the scenario doing the exponential transformation for all ten vertices, when only vertices change, only edge changes, both vertices and edges change table  presents the estimated computing time. it indicates that the proposed netdifm indeed runs fast, and the computational time increases as the network become larger.table computing time of four methods with sample size and permutations netdifm vewdm yates’d network ten vertices and edges; network vertices and edges; network vertices and edges the results of application network difference analysis for both the gwas of leprosy and gene expression data of ovarian cancer further confirm in practice that the proposed netdifm captured the network changes. shown in table  are the results of the proposed netdifm and other methods for detecting the network difference between two groups. the difference of gene interaction network with genes can be detected significantly at α =  by netdifm, netdifmpm and vewdm.table p-values of the four methods for the two real data sets netdifm vewdm yates’d group difference of the subnetwork of pi3k-akt signaling pathway was detected significantly at α =  by netdifm, netdifmpm and vewdm. when applied to the subnetwork of notch signaling pathway, all four methods can detect the network difference significantly. shown in table  are the results of the proposed netdifm and other methods for detecting the specific vertices, treating a vertex as well as its connected edge as a network, under permutation times.table p-values of the four methods for detecting the genes in the pi3k-akt signaling pathway and notch signaling pathway netdifm vewdm yates’d   kras   pik3ca   irs1   pdpk1   akt1   mdm2   tp53   dll3   dtx2   crebbp   ptcra   jag1   dvl2   snw1   hes1   rbpj   notch2   psenen   adam17   numb   ncor2 discussion and  complex disease is largely determined by a number of biomolecules interwoven into networks, rather than a single biomolecule. group-level comparison of network properties may shed light on underlying biological processes or disease mechanisms, and benefit the design of drug targets and drug combination for the therapy of complex diseases. meanwhile, although the conventional single-based paradigm has successfully identified a list of risk factors, one common sense is that there still exist an intermediate “black box” between the exposures and the disease phenotypes. in the “black box”, various risk factors weaved into complicated biological networks dominating the disease occurrence, development and prognosis. recent advances in high-throughput technologies and omics resources are revolutionizing biomedical research, and allow a transition from the traditional paradigm for biological and epidemiological studies of complex diseases to a new paradigm based on systems epidemiology. under this framework, network-based analysis has been integrated into observational study designs to organize the interdependencies of biomolecules and other factors at a human population level, expecting to open the “black box”. a key but inadequately addressed issue is still to develop valid statistical method to test possible differences of the networks between two groups. in our previous study, we have developed a statistical method for detecting the pathway effect contributing to disease, mainly under the framework of systems epidemiology. however, this method is limited to the pathway with chain structure, and can only capture the edge changes while omitting the vertex changes. meanwhile, the nonparametric bootstrap method has to be used to obtain the significance. at present study, we proposed a score-based powerful statistical test to detect the significant changes in biological networks between two different conditions. it can simultaneously capture the vertex changes and edge changes. various simulations were conducted to assess the reliability and statistical power of the proposed method. it indicated that both netdifm and vewdm were much more powerful than yates’d, and netdifm kept the best performance under various scenarios, and it can indeed capture the perturbation of vertices and edges in the network simultaneously. one strength for netdifm is that we can obtain its theoretic property, and thus can avoid the high computation burden. as expected, the proposed netdifm indeed runs fast. the vewdm was used fisher r-to-z transformation to identify significant differences between two correlations. fukushima et al. also developed an r package to identify differential correlations between two conditions based on fisher’s z-test which affords users a simple and effective framework in omics data. the vewdm can be treated as a global measure to detect the group difference of networks between two conditions, accounting for not only edges difference but also vertices difference. even though one is interested in testing particular vertex or edge rather than the whole network, its connected edge should also be considered. two real data sets analyses further highlighted that netdifm had more advantage in practice. in the gwas data of leprosy, we detected a candidate gene interaction network containing five genes. for the gene expression data of ovarian cancer, two candidate subnetworks, pi3k-akt signaling pathway and notch signaling pathway, respectively were considered and identified, suggesting that the proposed method is capable of identifying differential gene expression and gene-gene co-expression patterns, which are certainly helpful for us to further understand the underlying disease mechanism. rao et al. reported that combined overexpression of ova and mdm promotes oncogenesis by enhancing activation of the igf-1r–erk1/ signaling pathway, and jag enhances ovarian cancer cell growth and cisplatin-resistance. the expression of hes is confirmed to be strongly associated with the pathogenesis of ovarian endometriomas. meanwhile, decreased notch expression is associated with the poorly differentiated serous epithelial ovarian carcinoma histology. rbpj underexpression in ovarian tumor tissue relative to matched normal tissue. moreover, adam is one of the several metalloproteinases that play a key role in epidermal growth factor receptor signalling and can be a potential target antigen to devise novel immunotherapeutic strategies against ovarian cancer. the pi3k-akt and notch pathways and their abundant associated genes comprise complicated networks, which play a significant role in the progressive growth of tumor cells. network difference can result from not only changes of vertices but also changes of edges, and the changes of vertices-wise and edges-wise are often closely related. for instance, differential expression of genes may be due to either mutation of its own gene or the effects of expression changes of other genes in the network. however, the degree of differential expression of one gene due to its own mutation is often lower than affected by expressions of upstream genes in the network. reverter et al. presented an analytical procedure to simultaneously identify differential gene expression and connectivity for unweighted gene network. in their work, an edge between two genes is established if the absolute value of the correlation coefficient exceeds a fixed threshold. consequently, if we set the threshold less than, and the correlation coefficient between gene a and gene b is in cases and in controls, then the connection between gene a and gene b is treated as no difference between cases and controls. while in this situation, there exists a difference of the strength of connection between gene a and gene b among cases and controls, given our methods focus mainly on weighted biological networks. furthermore, the covariance structure between vertex changes and edge changes has been embedded into the proposed score-based network difference measure. in addition, one would be more interested in testing particular vertex or edge rather than the whole network or pathway. actually, a vertex as well as its connected edge can be treated as a subnetwork, and the proposed network difference method can easily be extended to identify the specific vertices. even though some local interventions were often generated to prevent and cure a particular disease, it is essential to understand the global system. the ‘think globally, act locally’ paradigm should be strongly embedded into our mind. one limitation in our paper is that we assume the network topological structure is fixed, and little attention has been paid on the network structure learning problem. constructing network structure means determining every possible edge with highest degree of data matching, and often one joint probability distribution of a number of variables can reflect more than one network structure. usually, combining experimental evidence with their experience, most biologists and clinical researchers have a growing awareness of the interplay between the biological components and can depict more or less the specific network or pathway for the corresponding biological process. meanwhile, numerous databases can be further borrowed to establish the network structure. the proposed netdifm will do not work in its current version under the scenario when the covariance matrix is not invertible. one possible solution is to first apply a shrinkage strategy to simplify the network, and then adopt the proposed statistic. for instance, we can first remove those edges if the correlation between the two linked vertices is smaller than a predefined threshold, and then apply the proposed test to the remaining network. statistical comparisons of group difference in biological networks are highly desirable. the proposed network difference measure was valid and powerful to detect biological network difference. r code to perform netdifm, netdifmpm and vewdm is provided in the additional file availability of supporting data the gene expression data of ovarian cancer were downloaded from the geo datasets. additional files additional file 1: table s the location and snp number for five susceptibility genes belonging to the network associated with leprosy. table s the location and probe sets number for genes belonging to the subnetwork from pi3k-akt signaling pathway and notch signaling pathway. table s type i error rates of four methods. table s type i error rates of four methods given the deviation from the normal distribution. figure s a network including vertices and edges. figure s the statistical power of the four methods under three scenarios. only vertices change, only edges change, both vertices and edges change. additional file 2: r code. abbreviations dcdifferentially connected dedifferentially expressed gwasgenome-wide association study keggkyoto encyclopedia of genes and genomes netdifmnetwork difference measure netdifmpmnetdifm based on permutation vewdmvertices and edges wise difference measure competing interests the authors declare that they have no competing interests. authors’ contributions jj, zy and fx conceived, designed the study and implemented the data analysis, jj, zy and xz drafted the manuscript. all authors read and approved the final manuscript. 
for additional gene expression data of ovarian cancer, two candidate subnetworks, pi3k-akt and notch signaling pathways, were considered and identified respectively. we therefore proposed a powerful score-based statistic to detect group difference in weighted networks, which simultaneously capture the vertex changes and edge changes. the proposed method, accounting for the vertex changes and edge changes simultaneously, is valid and powerful to capture the group difference of biological networks. group-level comparison of network properties may shed light on underlying disease mechanisms and benefit the design of drug targets for complex diseases. simulation studies indicated that the proposed network difference measure was stable and outperformed other methods existed, under various sample sizes and network topology structure. complex disease is largely determined by a number of biomolecules interwoven into networks, rather than a single biomolecule. a key but inadequately addressed issue is how to test possible differences of the networks between two groups. one application to real data about gwas of leprosy successfully identified the specific gene interaction network contributing to leprosy.
 complex disease is largely determined by a number of biomolecules interwoven into networks, rather than a single biomolecule. a key but inadequately addressed issue is how to test possible differences of the networks between two groups. group-level comparison of network properties may shed light on underlying disease mechanisms and benefit the design of drug targets for complex diseases. we therefore proposed a powerful score-based statistic to detect group difference in weighted networks, which simultaneously capture the vertex changes and edge changes.  simulation studies indicated that the proposed network difference measure was stable and outperformed other methods existed, under various sample sizes and network topology structure. one application to real data about gwas of leprosy successfully identified the specific gene interaction network contributing to leprosy. for additional gene expression data of ovarian cancer, two candidate subnetworks, pi3k-akt and notch signaling pathways, were considered and identified respectively.  the proposed method, accounting for the vertex changes and edge changes simultaneously, is valid and powerful to capture the group difference of biological networks. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords network medicinesystems epidemiologyscore-based statistical testnetwork comparisonhttp://dx.doi.org/ /501100001809national natural science foundation of china8157325931200994yuan zhongshang xue fuzhong issue-copyright-statement© the author 2016 