BACKGROUND
being a computational approach based on the prior knowledge, pathway enrichment analysis is widely used in the study of genotype-phenotype associations  <cit> .
biological pathway as a set of interactive genes  produces particular cellular response/outcome by executing a series of functional cascades.
it is curated by experts from wide range of science fields  <cit>  so that can supply more creditable functional details than general go module or network module.
different from exploring the unknown or indeterminate functions by network module, pathway-centered analysis always makes an effort to capture the permutation of established functions  in the change of phenotypes .
as a key approach of pathway-centered analysis, the pathway enrichment analysis or well-known gene set enrichment analysis   <cit>  can identify dysregulated pathway by qualitatively measuring the changed status of a pathway  <cit> .
in the pathway enrichment analysis, the dysregulation of a pathway is the most important issue  <cit> , and should be mathematically defined and measured well  <cit> .
it can estimate the conditional enrichment or status of a pathway, which is assumed to be associated with particular phenotypes.
current researches generally use genes with significantly differential expressions or differential correlations to evaluate the extent of the dysregulation of a pathway.
one kind of conventional method is evaluating the dysfunction of pathways in different conditions , such as fidepa   <cit> , spia   <cit>  and ipeap   <cit> .
the other kind is using pathways to characterize individual samples  <cit> , like corgs  <cit>  and pathifier  <cit> .
generally, all these methods focus on the genes with differential expression and their enrichments in pathways   <cit> , which assume the samples are of good purity in genotype-phenotype association study.
however, in the study of complicated phenotypes, e.g., cancer study, a relevant problem is the samples with the same disease phenotype might be full of different unknown subtypes due to disease heterogeneity  <cit> .
it is necessary to detect genes with new features observable in the complicated disease samples, and enhance the pathway enrichment analysis to be applicable in such previously unexpected situation  <cit> .
actually, there are new expression features extracted in recent studies, e.g., genes with differential expression variances  <cit> .
in the context of differential expression variance, it is still a challenge to measure the enrichment or status of a pathway.
a solution to this problem can promote the efficiency of pathway enrichment analysis on genotype-phenotype association because it will consider more complete information about the expression changes of pathway genes.
it can also provide new insights on the biological pathways by integrating additional expression and network features.
in this work, we propose a multiple-label based enrichment analysis to detect such dysregulated pathways, which simultaneously takes into account the genes with differential expression  and genes with differential expression variance  together .fig.
<dig> major differences between the measurements of dysregulated pathways used in conventional enrichment analysis and integrative enrichment analysis
obviously, the hypothesis underlying iea is that the dysregulated pathways involved in disease heterogeneity would be full of degs and/or devgs.
that means the identified pathways by iea would be disease pathways or their up-streams/down-streams .
however, current methods in pathway enrichment analysis only expect to give high-rank to disease pathways .
when iea identifies up-streams/down-streams of disease pathways, it further assistantly supplies a network of pathways to recover a global functional map and infer the associations among disease pathways and subtype-relevant pathways.
noted, the biological meaning of the edge in such network of pathways is the pathway crosstalk, which is just an important biological mechanism or functional relationship among pathways .
conventional researches tend to simply determine a pathway crosstalk by the overlapped genes in two pathways  <cit> , which disregard the statistical significance of the genes and interactions involved in the pathway crosstalk.
by contrast, degs and devgs in one pathway can be used as seeds, and further detected their interactive genes in the candidate crosstalking pathways by a random walk restart algorithm  <cit> .
the significance of a pathway crosstalk can be finally evaluated by the genes involved in this crosstalk as their enrichments in two pathways .
based on the above concepts and mathematical models, a new pathway-centered analysis framework, the integrative enrichment analysis , is implemented as  pathway enrichment score calculated by the hypergeometric test on differential genes ;  pathway crosstalk ranked by the random walk and hypergeometric test on rewired molecule networks;  pathway-phenotype association and subtype-factors determined by devgs in pathways.
according to a previously established evaluation scheme  <cit> , we found that, in particular cohorts , a few target disease pathways can be significantly high-ranked by iea, which supplied the evidences of the deviation-based disease characteristics , and iea is more effective than other state-of-the-art methods in this condition.
furthermore, by a proof-of-concept study, we shows the details of iea on analyzing real transcriptional data related to complex diseases, e.g., diabetes and colorectal cancer.
iea indeed captures the previously under-estimated pathways full of devgs and degs.
these newly identified dysregulated pathways would be heterogeneity-relevant pathways and are found to be significantly linked to disease pathways  by estimated crosstalks.
many candidate subtype-factors are also recognized as devgs or pathways associated with the risk of subtypes of genotype-phenotype associations.
totally, iea supplies a new way of over-representation approach  <cit>  to carry on enrichment analysis in the complicate context of clinical application , and could be easily expanded to functional class scoring or pathway topology based approaches , which will be a necessary complementary and cooperative approach to conventional ones  <cit> .
the matlab scripts of the software named ieapackage and some alternative r scripts have been deposited in github and accessed in https://github.com/bluesky2009/integrative-enrichment-analysis.
this software has been developed and tested in windows  <dig> or windows  <dig>  and matlab  <dig> or matlab  <dig>
methods
generally, enrichment analysis includes three categories of methods: over-representation approach, functional class scoring and pathway topology based approaches.
although these methods are all focusing on evaluating the phenotype-associated pathway, they would be based on different hypothesis.
this work and the proof-of-concept study are based on the over-representation approach, which measures the dysregulation extent of a pathway according to the number of dysregulated genes in this pathway.
traditional methods only evaluated the degs in a pathway; by contrast, iea evaluates the degs and devgs in a pathway.
thus, the meaning of the statistic for the integration of iea is as completely as possible to measure the dysregulation extent of a pathway according to the number of dysregulated genes  in this pathway, which have been well defined and introduced in follows.
differential gene expression and differential expression variance
given a gene x has expression profiles in control and case samples as x and x’ respectively, the expression variance of this gene in control and case condition are e2) and e2) respectively.
here, u and u’ are average expressions of gene x in control and case samples respectively.
then, the conventional criterion and measurement of genes with differential expression  are:  <dig> h0:ex=ex’;h0rejected; where x or x’ are the original/raw expression levels.
noted, the differential expression includes up-regulation  and down-regulation .
except for these degs , the genes with differential expression variance are also discriminative features  <cit> .
the expression variance concerned features, e.g., bimodal gene expression, is already known as an important expression pattern in the control of a transition of biological systems  <cit> , such as: disease development, cellular differentiation, and phase transition.
however, the differential expression variance of genes has not been studied in a systematic way to the best of our knowledge, especially for its usage in the pathway enrichment analysis.
the differential expression of genes, used in conventional enrichment analysis, requires the gene’s expressions under different conditions to distribute around different mean expression levels .
by contrast, differential expression variance of genes  can be defined as the genes’ deviations being significantly different under dissimilar conditions , such as:  <dig> h0:ex‐u=ex’‐u’;h0rejected;andh0:ex=ex’;h0notrejected where x or x’ is the original expression level, |x-u| or |x’-u’| is the relative expression level.
noted, the differential expression variance includes tight-regulation  and relax-regulation .
and importantly, as defined above, the devgs have excluded degs, or there is no overlap between devgs and degs in this work.
that means, when one gene has both differential expression and differential expression variance, this gene is thought as deg in priority in order to be consistent with conventional analysis; and, of course, this kind of genes are worthy of deep research in future work.
actually, given x or x’ satisfy normal distribution, |x-u| or |x’-u’| will be folded normal distribution, then the wilcoxon rank sum test instead of student’s t-test is used in the significance test of devgs.
integrative enrichment analysis in the context of differential expression variance
obviously, the conventional enrichment analysis limits to estimate the extent of differential expression rather than differential expression variance.
when considering the contribution of devgs on pathway’s dysregulation, it is necessary to refine the conventional approach to take into account the degs and devgs together.
naturally, an easiest strategy is to put degs and devgs together as the same dysregulated genes and use conventional hypergeometric test to obtain the p-value.
however, this will disregard the respective distribution of degs and devgs in a target pathway and in the whole transcriptome.
thus, we extended the hypergeometric test on two kinds of enriched genes simultaneously as bellows.
our approach, noted as ht <dig> , still depends on the hypergeometric distribution and uses p-value to measure the dysregulation of a pathway in the context of differential expression variance.
briefly seen in table  <dig>  given there are expression data on total n genes, and x <dig> degs and x <dig> devgs selected respectively.
for some pathway, k <dig> and k <dig> genes from pathway members  have differential expression and differential expression variance respectively.
then the significance of deregulated genes as degs or devgs enriched in this pathway can be estimated by formula  <dig>  this p-value also ranges from zero to one.
the less the p-value is, the larger dysregulation extent the pathway has, when the significantly larger number of genes in this pathway show differential expression or differential expression variance.
<dig> px1=k <dig> x2=k2=x1k1x2k2n−x1−x2y−k1−k2nypx1>k <dig> x2>k2=1−∑<i <dig> i2>∈0x1×0x2−k <dig> x1× and disease-relevant pathways .
given several genes in a pathway as seeds, iea uses random walk to find their partner genes in the other pathway.
in fact, random walk with restart  is a well-known ranking algorithm for candidate gene prioritization  <cit> .
it supplies the probability of searching the random walker at nodes in the steady state, so that, it can give a measure of proximity between source nodes  and other nodes in molecule network .
let n be the adjacency matrix of a gene network with node set v and edge set e, in which the element nij equals one if e ∈ e  represents the interaction between genes/nodes i and j), or zero otherwise.
based on the topological structure of the gene network, the transition matrix t can be calculated.
each element in the transition matrix is denoted as tij and represents the probability of transition from node i to node j. the value of tij can be given by one of two ways as follows, the first one is topology-weighted and the second one is correlation-weighted.
tij=nijdi,ifeij∈e <dig> otherwise,wheredi=∑j∈vnij  tij=wijnijwi,ifeij∈e <dig> otherwise,wherewi=∑j∈vwijnij the rwr algorithm  <cit>  updates the probability vectors by pk+1=1−λtpk+λp <dig> k> <dig> where t is the transition matrix and p <dig> is the initial probability vector with the sum of the probabilities as one.
in p <dig>  all the source nodes are assigned equal probabilities and other nodes are given zero.
p∞ is obtained when the algorithm is convergent.
if p∞ > p∞, node i is thought to be more proximate to source nodes than node j does.
thus, a two-way rwr approach  is proposed to search the genes involved in two interactive pathways and estimate their enrichment for evaluating the pathway crosstalk.
the steps of two-way rwr include: for each pathway u, its degs and devgs are used as source nodes/genes, and rwr is used to rank the genes in known molecule network, e.g., protein association network collected from string database  <cit> .
in the high-ranked genes from above rwr analysis, the genes belonging to pathway v are the partner genes interactive with source genes.
based on the sources genes and their partner genes, the enrichment of those interactive genes  in pathways u and v can be evaluated by our ht <dig> approach, i.e., p-value in formula  <dig>
for every pathway, the analysis in steps  and  is repeated.
then, given a pathway pair , it is a pathway crosstalk only when euv and evu are both significant.
finally, the map of pathways consist of those selected pathway crosstalks, where a node represents a pathway and an edge represents a pathway crosstalk.
screening subtype-factor of genotype-phenotype associations based on devgs and dysregulated pathways supplied by iea
the second assistant down-stream analysis method of iea is to screen subtype-factors according to the available clinical indices.
as stated above, iea focus on the devgs and their involved pathways, and these genes and pathways are thought as signatures of potential subtypes of heterogeneous samples.
however, these hidden subtypes might have not been identified or formalized in clinics.
to evaluate such new signatures or subtypes, one direct strategy is to measure the correlation between genetic signatures  and clinical indices .
if one signature is significantly related to some clinical index, the subtype represented by such signature would be medical meaningful as to be observable in clinics and this signature is also called as subtype-factor related to particular clinical index.
the approach to identify such subtype-factors is described in bellows.
for each pathway, its devgs are used to group case  samples into two clusters, when the case  samples have high varying expression compared to control  samples.
that means these genes have over-expression in one group of samples and under-expression in the other group of samples.
this pathway would be a candidate subtype-factor when these two sample clusters are discriminative on some clinical index.
on this condition, a clinical subtype of samples is thought to be related to a given clinical index, which is represented by a subtype-factor .
obviously, the clinical subtype of a particular sample might be contributed by many subtype-factors .
given a known phenotype , a few subtype-factors correlated with this phenotype can be found, although which just reveals only the tip of the iceberg for the subtypes of genotype-phenotype associations.
particularly, different from conventional un-supervised clustering for subtype identification, a supervised-like clustering approach  is proposed to identify subtype-factors on the level of pathways.
firstly, the case samples can be grouped into two clusters according to their features’ values  compared to those values of control samples: on each feature , one group of samples have larger values than controls meanwhile the other group of samples have less values than the same controls, or vice versa.
that means, a hyperplane determined by a few control samples could separate the samples space into two sub-spaces, and case samples in each of two sub-spaces are grouped into one cluster.
secondly, some clinical information of samples can be used to evaluate the potential subtype represented by such two clusters of case samples.
if the clinical values of these two groups of samples have significant difference, a clinical subtype of genotype-phenotype association  is identified and the corresponding pathway is a subtype-factor corresponding to the given clinical index.
practically, the slc algorithm on a pathway is implemented as bellows: discrete the expressions of devgs of case samples into binary vector based on the values of controls: for a devg, if its expression value is larger than the mean of controls, it is one in the binary vector; otherwise, it is zero.
clustering case samples based on the binary vectors by conventional methods as hierarchical clustering or k-means, which obtains two sample clusters.
calculating the significance of difference between clinical indices among above two sample clusters.
if the difference is significant, this pathway is identified as a subtype-factor of the association between the given pathway and clinical index.
RESULTS
the evaluation of biological meaning of iea by method comparison
iea is proposed to evaluate dysregulated pathways by differential gene expression and differential expression variance together.
differential expression variance has been reported as a new and important expression change during a phenotype change  <cit> , e.g., diseases.
in this work, the biological hypotheses underlying iea is that, the dysregulated pathways full of genes with differential expression variance would be subtype-relevant pathways.
although subtype-relevant pathways for particular complex disease are unclear in current pathway databases, e.g., kegg, it is still able to investigate if prior-known disease pathways in kegg would be subtype-relevant and if iea can identify them.
in the previous study of gene-set analysis  <cit> , a comparison scheme has been built to evaluate the performances of different enrichment analysis methods  based on multiple expression datasets about complex diseases.
different from previous general comparison, we focus on the comparisons by approach-specific datasets, in order to mainly evaluate the biological meaning of iea.
according to the comparison protocol  <cit> , we ran total eight representative enrichment analysis methods on  <dig> geo datasets with target pathways in kegg, and obtained the rank of target pathway estimated by each method on each dataset; then, for each dataset, we rank the eight methods according to their prioritization performance or sensitivity performance  <cit> , and this dataset is assigned as a specific-data for the top-k methods ; thus, all specific-data for one method can consist of k-order approach-specific dataset.
generally, on one method’s approach-specific dataset, this method should have best or comparable performances than other methods, so that, the biological characteristics assumed by this given method would significantly displayed on these datasets.
therefore, we can use this strategy to investigate the biological meaning of iea in real datasets.
in bellows, we firstly summarize the biological hypothesis hold by different state-of-the-art enrichment analysis methods and their respective quantitative measurements, and then discuss the comparison between iea and others.
plage: it assumes the activity of pathway rather than the expression of pathway genes determines the activated or inhibited status of pathways under different conditions; and the pathway activity is measured by an activity score as the weights of a metagene extracted from all pathway genes by svd   <cit> .
gsva: it proposes the change of pathway activity between control and case should be evaluated at the level of samples, e.g., considering the variation of pathway activity over a sample population; and the pathway activity is measured by so-called gsva score as a function of the expressions of genes inside and outside the pathway, and these scores are assessed similarly as gsea by using the kolmogorov-smirnov  like random walk statistic  <cit> .
padog: it assumes that, if the genes highly specific to a given pathway occur differential expressions, the respective pathway would be truly relevant in that condition; thus, a new gene set score is calculated as the mean of absolute values of weighted moderated gene t-scores where the gene weights are designed to be large for the genes appearing in few pathways and small for genes that appear in many pathways  <cit> .
globaltest: it holds an assumption that, if a group of genes  can be used to predict the clinical outcome, the expression patterns of such gene group must differ for dissimilar clinical outcomes; thus, it uses generalized linear model to give one p-value for a group of genes, not a p-value for each gene, which can be applied to estimate the enrichment of a given pathway  <cit> .
mrgse: it proposes that the high ranks of expression changes  of genes can indicate the differential expression of a set of genes ; and the enrichment score or the test statistic of a pathway is the mean rank of this gene set, i.e., the average of the ranks of t-statistics of pathway genes  <cit> .
gsa: it is similar to gsea, and proposes two improvements as the maximal average statistic for summarizing gene-sets, and restandardization for accurate enrichment inferences  <cit> .
ora: it takes into account the number of differentially expressed genes observed in a pathway as indicators of pathway states; generally, it uses a basic contingency table to test the association between the differential expression status of a gene  and its membership in a given gene set , which can be measured by the p-value of a hypergeometric test  <cit> .
iea: it is proposed in this work to generally consider the contribution of expression variance in a dysregulated pathway; as one implementation, this work takes into account the number of degs and devgs observed in a pathway as indicators of pathway states; it is designed to test the association between the differential expression/differential expression variance status of a gene and their memberships in a given gene set, which can be measured by the p-value from proposed ht <dig> approach in this work.
first of all, we can cluster the above eight approaches by their performances on all datasets to investigate the general association among different methods.
as shown in figs.
<dig> and  <dig>  the similarity among any two methods is measured by four kinds of criterion: the first one is whether the ranks given by two methods on the same dataset are also the same ; the second one is whether the ranks given by two methods have the same change tendency among different datasets ; the third one is whether the p-values given by two methods on the same dataset are also the same ; and the last one is whether the p-values given by two methods have the same change tendency among different datasets .
obviously, gsa and padog are both based on conventional gsea, so that they are similar; the proposed iea is based on ora, thus, they also have similar performances on different datasets; plage and globaltest are closely clustered together, one reason is that they both estimate a score from all pathway genes rather than individual genes ; in addition, mrgse and gsva are much different, and also different form other methods, which is possibly because they have specific design principles on the measurement of pathway dysfunctions, i.e., mrgse combines the t-statistics of individual pathway genes meanwhile gsva uses a score as a function of the expressions of genes inside and outside a pathway.fig.
<dig> category of representative gene set analysis approaches based on clustering of prioritization performance.
a method clustering based on euclidean distance of ranks of all pathways.
b method clustering based on correlation distance of ranks of all pathways
fig.
<dig> category of representative gene set analysis approaches based on clustering of sensitivity performance.
a method clustering based on euclidean distance of p-values of all pathways.
b method clustering based on correlation distance of p-values of all pathways
then, we directly grouped the datasets according to the performance of a given method, e.g., some datasets are included as k-order iea-specific datasets, only when the rank of iea performance compared to all methods are in the top-k on these datasets, where k is set  <dig> in this study.
to quantify the performance, sensitivity  and prioritization  are adopted as previously  <cit> .
in previous evaluation on these datasets, padog displays consistently comparable performance with other methods, meanwhile, plage, globaltest and mrgse have the best performances on some categorise of datasets  <cit> , which already suggest the existence of approach preferences.
thus, the categories of datasets induced in this work can indicate potential preferences of different methods on particular datasets.
of course, on the preference-specific datasets , the preferred method should have the best performance; meanwhile, some other methods would have comparable performances.
therefore, this new evaluation scheme can supply evidences for two hypotheses: one is that the expression features of a pathway measured by a given method actually have biological meaning, and are existing or observable in real datasets corresponding to particular phenotypes ; the other one is that a method has significant preferences  on a group of datasets , which is comparable or not by other methods.
as seen in table  <dig> for prioritization performance, every method shows its preferences on particular datasets ; in addition, the method comparison on more strict specificity as 2-order approach-specific datasets  and more weak specificity as 4-order approach-specific datasets  have also been done and reported in additional file 2: table s <dig> and additional file 3: table s <dig>  combined these results together, it is easily to see that: iea, ora and padog have significant preferences, due to their highest performances than all other methods on their correspondingly preferred datasets; globaltest, gsva, plage and mrgse also have significant preferences on datasets, although sometimes padog could have comparable performances with them ; besides, gsa is challenged by padog, because on the gsa-specific datasets, padog always achieves similar or better performances, but not vice versa.table  <dig> the prioritization performance about method comparison on approach-specific datasets
gsa
padog
iea
mrgse
ora
globaltest
gsva
plage
the performance of an approach on its specific dataset is highlighted in bold.
and the performance of comparable approaches on some specific dataset is highlighted in bolditalic
and seen in table  <dig> and additional file 4: table s <dig> for sensitivity performance, again, every method shows its preference on particular datasets .
it seems that globaltest and plage have generally comparable performances with other methods according to their performances on many preferred datasets of other methods.
even though, iea shows the best performance on the iea-specific datasets, and is better than ora on the ora-specific datasets.
this fact strongly suggests that, iea actually can detect dysregulated pathways, and displays competitive performance than conventional methods in current evaluation scheme when those dysregulated pathways are just target pathways ; besides, iea is realized based on the conventional ora, and improves ora on the sensitivity performance, which would just be contributed by considering the new feature genes as devgs in dysregulated pathways.table  <dig> the sensitivity performance about method comparison on approach-specific datasets
gsa
padog
iea
mrgse
ora
globaltest
gsva
plage
the performance of an approach on its specific dataset is highlighted in bold.
and the performance of comparable approaches on some specific dataset is highlighted in bolditalic
finally, above dataset-driven method comparison supplies new insights on the performance specificities of many gene-set approaches, especially for iea; and also supports the importance and biological meaning of dysregulated pathways identified by iea.
we can draw following conclusions: although a few methods have consistent performances on many datasets , different method still have their preferences on the expression characteristics of dysregulated pathways, so that each method can achieve significantly better performances on their specific datasets rather than all datasets.
especially, on the iea-specific datasets, iea indeed are the best one than other methods; and even on the ora-specific datasets, iea is better than ora on sensitivity performance.
on the iea-specific datasets, the target pathways or disease pathways are possibly just the subtype-relevant pathways, so that, iea have competitive performance in the comparison scheme.
therefore, the expression variance focused by iea is actually full of biological meaning, and will help iea to detect new dysregulated pathways, e.g., subtype-relevant pathways.
in addition, the complex diseases concerned in iea-specific datasets actually already have reports about the existence of subtypes on genetic level, such as dilated cardiomyopathy  <cit> , renal cancer  <cit> , prostate cancer  <cit> , colorectal cancer  <cit> , and thyroid cancer  <cit> .
every method, or every method category, can actually capture particular dysregulated pathways.
when the target pathway of a dataset just displays the expression characteristics focused by an approach, such approach would have better performance on this dataset.
on some specific datasets preferred by other methods, iea should face two conditions: one is under the condition that the target pathway is not subtype-relevant pathway, and iea will be underestimated but has supplied a useful down-stream analysis  to assistantly link the target pathway and potential subtype-relevant pathway identified; the other one is under the condition that the target pathway is a subtype-relevant pathway, and iea should be further enhanced by integrating expression variance with other pathway pattern , which is worthy of study in future.
a proof-of-concept study of iea on transcriptional analysis of complex diseases
iea has been applied to detect the biological malfunction of complex diseases  on the pathway level rather than gene level.
iea, as a pathway-centred analysis approach, not only supplies the conventional pathway enrichment analysis but also extracts divergent pathway associations, such as: pathway & disease genes, pathway & pathway , and subtypes of pathway & clinic .
firstly, it is the data pre-procession.
data needed in iea have been prepared from the public resources: the gene expression data of human islets from non-diabetic and diabetic were downloaded through geo  <cit> ; there are two datasets, the main dataset gse <dig>  <cit>  contains samples from  <dig> non-diabetic and  <dig> diabetic with  <dig> genes, and the replicate dataset gse <dig>  <cit>  contains samples from  <dig> non-diabetic and  <dig> diabetic with  <dig> genes; the gene lists of  <dig> kegg pathways are obtained from gsea package  <cit> ; the human protein interaction network  are extracted from string database  <cit>  with confidence score no less than  <dig> ; four clinical indices are also obtained from the supplementary of original study  <cit>  as sex, age, bmi, and hba1c; meanwhile, diabetes associated genes are searched from genecards database  <cit> .
secondly, it is the main step of iea.
different scores of pathway enrichment are calculated:  the conventional score  as p-values of hypergeometric distribution of degs in a pathway;  the conventional score  as p-values of estimated pathway enrichment;  the new score  as p-values of hypergeometric distribution of differential genes  in a pathway calculated by the proposed ht <dig> approach.
for ora or iea, the thresholds of p-value of significance test on degs or devgs are both set as  <dig> , and adopted a pervious strategy to select those feature genes  <cit> : 1) select all genes with fdr adjusted p-values no more than  <dig> ; 2) if the genes selected are less than  <dig>  re-select all genes with p-values no more than  <dig>  and fold-change no less than  <dig> ; 3) if the genes selected are still less than  <dig>  directly use the top 1 % of genes ranked by p-values from least to largest.
thirdly, it is one assistant step of iea.
pathway crosstalks are evaluated by two-way rwr approach.
the interactions selected from pin consist of differential network  <cit> , where the selected interactions have significant correlation difference between diabetic and non-diabetic groups.
on this differential network, two-way rwr approach is used to find the pathway crosstalks.
the most significant pathway crosstalks  consist of the map of pathways.
besides, the enrichments of pathway genes or disease-associated genes in the high-ranked genes of rwr are also analyzed and evaluated by auc  <cit> , which support the efficiency of rwr on pathway-related analysis.
fourthly, it is the other assistant step of iea.
the devgs in each pathway are used to group samples in two clusters by slc approach.
for each clinical index from sex, age, bmi, and hba1c, its subtype-factors are identified.
as comparisons, the all genes of each pathway are also used to directly group samples to evaluate the significance of detected subtype-factors.
obviously, the above analysis routine can be applied on any other dataset of samples with different kinds of complex diseases.
noted, our proof-of-concept study is to combine a group of genes with differential expression and a group of genes with differential expression variance.
we don’t select the genes with high variances, but the genes whose expression variances can distinguish different conditions/phenotypes.
if the expression variance of a gene can classify samples with different phenotypes well, the dominant component of this gene’s variance could be biological variance.
obviously, the technique variance should have no such discrimination.
we have also checked the correlation between the variance of each devg and each clinical index.
many genes’ expression variances even have significant correlation with clinics, which would not be caused by technique variance too.
our experiment on microarray is an application of iea to support the idea of combining deg and devg, and the results show our method’s efficiency.
indeed, many approaches are still proposing to improve the selection of conventional deg or even new devg.
the removal of technique variance will improve the selection of devg and final iea, which is our future work.
diabetes associated genes on pathways
first of all, we investigated the overlap between prior-known disease genes  and pathway genes.
many pathways are full of diabetes associated genes , which means pathways could have great changes during disease development and progression.
these pathways would be causes or outcomes of the disease.
although iea pays attention to the identification of dysfunctional pathways  rather than discrimination of causal pathways, as introduced in follows, the map of pathways can further complementally supply some clues of the causal roles of pathways at the level of network of networks  <cit> .
then on the main dataset gse <dig>  <cit> , we have obtained feature genes as summarized in table  <dig>  there are  <dig> degs and  <dig> devgs selected by iea, many of them are also detected on the replicate dataset gse <dig>  <cit> .
the  <dig> genes of degs are disease genes, and  <dig> genes are pathway genes; meanwhile, the  <dig> genes of devgs are disease genes and  <dig> genes are pathway genes.
obviously, there are many disease-informative or function-informative genes disregarded in conventional analysis , and iea can capture these genes and estimate their effects in the dysfunction of pathways.table  <dig> the statistic on degs, devgs and their overlapping with pathway or disease genes
adeg points genes with differential expression; devg points genes with differential expression variance; pg_deg points the pathway genes in degs, i.e., the overlaps between pathway genes and degs; dg_deg points the disease genes in degs, i.e., the overlaps between disease genes and degs; pg_devg points the pathway genes in devgs, i.e., the overlaps between pathway genes and devgs; dg_devg points the disease genes in devgs, i.e., the overlaps between disease genes and devgs; degup and degdown point genes with up-regulation and down-regulation respectively; devgup and devgdown point genes with relax-regulation and tight-regulation respectively
furthermore, in degs, there are  <dig> gene up-regulated in disease state and  <dig> genes down-regulated.
meanwhile, there are  <dig> genes tight-regulated in disease condition and  <dig> genes relax-regulated.
the examples of such four expression patterns are shown in fig.
<dig>  myc, known as a cancer oncogene, is also reported to be altered in diabetes  <cit> .
seeing fig.
4a, myc is indeed up-regulated in the diabetes state.
insulin is known as a main cause of diabetes  <cit> , and its an isoform as ins-igf <dig> actually has down-regulation when diabetes occurs .
hoxd <dig>  as a gene in the homeobox family encoding a highly conserved family of transcription factors, has an important role in the morphogenesis.
it has tight-regulation in diabetes condition , so that, it may participate in the accurately regulation  <cit>  of biological processes associated to diabetes.
by contrast, rexo <dig>  known as transcription elongation factor b polypeptide 3-binding protein  <dig>  would be a cofactor involved in gene regulation  <cit> .
this gene shows relax-regulation in diabetes condition .
thus, rexo <dig> would be a cause or indicator of some subtypes of diabetes.
indeed, the original study has supplied four clinical indices  <cit> , we found age is mostly related to the subtypes or sample clusters of diabetes determined by rexo <dig> .fig.
<dig> the cases of differential expression patterns of deg and devg.
for genes myc, ins-igf <dig>  hoxd <dig>  and rexo <dig>  - give their expression profiles on dataset gse41762; - give their expression distribution on dataset gse41762; - give their expression profiles on dataset gse38642; - give their expression distribution on dataset gse38642
more importantly, devgs  don’t mean no-changes.
as the key point of our model and method, devg means a gene would denote activation of a signalling pathway  in a group of samples, meanwhile, inactivation of this signalling pathway  in another group of samples.
this would be a main cause of heterogeneous samples.
the biological mechanism underlying this phenomenon would be the switch of pathway activation.
if based on other kinds of enrichment analysis framework, it is possible to discuss the activation, inactivation, or activation-switch of a signalling pathway, which will be studied in our future work.
dysregulated pathways identified to capture degs and devgs simultaneously
in conventional analysis, the genes with differential expression are focused; now, the genes with differential expression variance are also attractive.
in the context of differential expression variance, the dysregulated pathways are expected to have as many degs & devgs as possible, which can be captured by over-representation approach like the proposed iea.
for evaluation, a gene-distribution graph is further introduced to show the percentages of degs and devgs respectively for each pathway.
in fig.
<dig>  a pathway is represented by a point whose indices in axis are the percentages of degs and devgs in this pathway respectively.
obviously, the pathways full of degs and devgs tend to locate at the right-up conner of such gene-distribution graph.fig.
<dig> the gene-distribution graph of pathways identified by different methods.
a-c the pathways high-ranked  by ora, iea and gsea respectively, whose results on dataset gse <dig>  d-f the pathways high-ranked by ora, iea and gsea respectively, whose results on dataset gse38642
to validate the pathways identified by iea and compare to other methods ora and gsea, we firstly calculated the percentages of degs and devgs of each kegg pathway and plot them on the gene-distribution graph.
then we calculated the dysfunctional score  and rank all kegg pathways, and found the top- <dig> selected pathway  in graph.
obviously, the pathways high-ranked by ora have least scores located at the right-bottom of gene-distribution graph .
it means the pathways selected by ora are full of degs rather than devgs.
by contrast, in the iea, the detected pathways are full of degs and devgs respectively and significantly .
even more, the dysregulated pathways detected by well-known gsea  <cit>  on these datasets shown weak performance on the identification of pathways full of devgs .
thus, iea indeed can effectively detect the pathways under-scored in conventional analysis, and these dysfunctional pathways would be disease-relevant or subtype-relevant.
this conclusion is also supported by the similar results from the analysis on replicate dataset .
noted, the pathways identified by different methods can be significantly observed in the analysis on the replicated dataset .
the issue of the consensus of pathway identification is not discussed more here, although some other studies have worked to improve the robustness of dysregulated pathway identification by integrating other prior information   <cit> .
in the high-ranked pathways identified by iea, many pathways are actually full of devgs.
these devgs could be further associated to potential subtypes of samples, which are discussed in follows.
besides, some of these pathways under-scored by other methods indeed have been reported to be altered in the disease state .
for examples, 'kegg hematopoietic cell lineage'.
diabetes is known to compromise the function of the bone marrow   <cit> , and diabetic complications mainly including macrovascular events might be from the dysfunctional bm-derived hematopoietic cells.
'kegg cytokine-cytokine receptor interaction'.
cytokines regulate inflammatory and immune responses, which play important roles in the pathogenesis of diabetes and its microvascular complications.
the functional variations of cytokines and their receptors can benefit the prediction of the susceptibility and progression to diabetic nephropathy   <cit> .
as the potential pathogenic mediators in dn, cytokines might provide new potential therapeutic agents for disease treatments.
dysregulated pathway-crosstalk identified to reveal the interactive map and module among pathways
different from conventional strategy to use the overlapping genes as pathway crosstalks  <cit> , two-way rwr is assistantly applied to find the interactive genes between any two pathways.
rwr is previously used in the rank of disease genes  <cit> , which holds the assumption that the candidate pathogen genes are more proximate to the known disease genes than randomly selected genes.
this assumption is expected to be hold for pathway genes too.
we first evaluated the pathway genes possibly selected by rwr,which uses the identified degs & devgs in a pathway as seeds.
in the high-scored genes by rwr, there is a significant amount of pathway genes , which is the same as disease genes ranked by rwr .
in the evaluation in table  <dig>  we have set two kinds of control experiments.
one is the prior-known network used, i.e., two sources as string  <cit>  and hprd  <cit>  are both applied; the other one is the feature genes used in the given background network, i.e., three kinds of feature genes  are respectively used to calculate the auc values  <cit>  to evaluate the efficiency of selecting/ranking pathway genes or disease-associated genes.
depending on these experiments, we can find:table  <dig> the auc of different rank lists for pathway genes
afor feature genes like all genes and degs, they don’t use network information, so that, they have the same auc values on the same dataset although different network used
afor feature genes like all genes and degs, they don’t use network information, so that, they have the same auc values on the same dataset although different network used
as a control, when all genes are ranked according to p-values, its auc is low.
and when the selected degs are ranked according to p-values, the auc is similar to that of all genes.
meanwhile, when the genes from two-way rwr excluding seeds are ranked according to proximity values, the auc achieves highest, which support again rwr-based approach is effective to capture interactive phenotypic genes as pathway genes or disease-associated genes.
thus, two-way rwr is effective to mimic the pathway crosstalk and construct the associations among pathways, which is obviously consistent in multiple control experiments by using different prior-known network , or different datasets , or even different ranked/selected feature genes .
the protein association network  rather than protein physical network  would be more efficient to lead the two-way rwr to link seed genes to pathway genes or disease-associated genes.
there are at least two reasons for this result: one is that the known protein physical network is still greatly incomplete, by contrast protein association network would supply additional predicted interactions with high confidence; the other one is that, except for direct interactions between pathways, protein association network would cover much more indirect interactions or long-term interactions, which would mimic the pathway crosstalk well.
in all, protein association network is efficient to detect the associations among pathways, whose usage to accurately predict physical interaction of pathways is out of this work and would be a future study.
thus, rwr is actually effective to detect the genes interactive within known pathway genes.
then by two-way rwr, we can find the interactive genes from two pathways, and select any pathway-pair as a crosstalk significantly .
all the crosstalks connect the known pathways as a map, where each crosstalk is an edge and a pathway is a node.
this map of pathways is a network of networks, rather than the original background network of separate molecules.
in the map of pathways, the modules of interactive pathways can be detected, where a module represents a group of closely inter-connected pathways.
in these pathway modules, two modules are obviously related to diabetes.
one module  is the group of signalling pathways.
signalling pathways are known as the up-stream functions in the cascades of signals, so that, they have great possibility to be the causes of the dysfunction of down-stream functions, e.g., diabetes pathways.
the other module  is just the group of pathways concerning glycolysis, sugar metabolism, glycosaminoglycan and diabetes, which seems to be extremely a core pathway module of diabetes.
noted, the type ii diabetes pathway would have significant crosstalk with pentose phosphate pathway as shown in fig.
<dig>  as reported, the pentose phosphate pathway is widely activated in diabetes and its complications , thus this pathway would be important to understand the risk of diabetes diagnosis and treatment in clinical application.
obviously, the pentose phosphate pathway have no significance on the selection by iea or other methods, and actually, its importance is reflected from the topological structure of the map of pathways.
these facts reveal:  conventional approaches usually focus on single pathways, so that, they can sometimes capture the disease associated pathways relevant to particular phenotypes but can’t distinguish or underestimate the potential causal relationship among pathways;  the proposed iea supplies the map of pathways to reflect the functional organization of pathways, and disclose the key modules of pathways, such as the upstream pathways related to signalling pathways and the downstream pathways associated to diseases;  on the map of pathways, those pathways full of devgs and degs tend to interact with disease pathways, indicating the determinant of subtypes  are actually also the potential determinants of diseases.
totally, the map of pathways supplies us a new viewpoint of functional organization at the level of network of networks.fig.
<dig> the main topological structure of the map of pathways reconstrcucted on diabetes datasets.
module  <dig> is full of signalling pathways, and module  <dig> is full of disease/diabetes associated pathways
dysregulated pathways associated to clinical indicator as subtype-factors
the pathways full of degs and devgs are specially selected by iea, in which the devgs might be the cause of potential subtypes of samples.
compared to the mean values of genes in case, the samples in control can be divided into two groups: one group of samples have larger gene expressions than that mean value; and the other group of samples have fewer values, or vice versa.
thus, we have applied two strategies to cluster the samples in control or case group, which can associate one pathway to some clinical index.
the first common strategy  is using the expressions of whole genes in a pathway to cluster samples into two clusters, and test the significance of these two groups of samples on one clinical indicator .
the second strategy  is using the discrete value of devgs in the same pathway to group samples and measuring their relationships with clinical indices.
additional file 9: table s <dig> gives the p value of the association of each pair between a pathway and a clinic index for normal samples, and additional file 10: table s <dig> gives those for diabetes.
generally, slc tends to discover more significant potential subtypes of samples corresponding to particular clinic index .
thus, devgs actually have more power to identify the subtypes of genotype-phenotype associations than conventional approaches based on differential expression only.fig.
<dig> the summary of dysregulated pathways identified as subtype-factors corresponding to four clinical indices.
for example, there are  <dig> pathways identified by iea to be correlated with hba values , but, there are only  <dig> pathways recognized by conventional strategy .
similarly, for other clinical indices, iea can also identify many subtype-factors
obviously, one pathway can associate to multiple clinical indices, and one clinical index can also relate to multiple pathways.
‘kegg wnt signaling pathway’ is consistently related to bmi in the analysis of duplicate datasets.
in fact, this pathway already has a few evidences on their associations with body weight .
‘kegg cell cycle’ is found to possibly associate with age, sex, bmi or hba, which is known as a common factor in disease development and progression  <cit> .
although each devg can separately determine some potential subtype-factor, here, we only discuss the determination of devg combination at the level of pathway.
in addition, the association of each pair between a gene  and a clinical index are listed in additional file 11: table s <dig> and additional file 12: table s <dig> for normal and diseased samples respectively.
noted, subtype is a potential biological explanations of devg, so that, we have evaluated the possible subtype-factors rather than subtypes by the association between devg  and known clinical index.
in the statistical analysis, we have analyzed age, sex, bmi, and even hba1c which evaluates the risk of diabetes.
thus, this clinical information can help us to find subtype of genotype-phenotype associations as stated above.
of course, in breast cancer, there are some well-known subtypes determined on genotypes.
to the best of our knowledge, in diabetes, the subtypes as t1d, t2d, gestational diabetes, surgically induced diabetes, chemically induced diabetes, are not defined by one or two genes/proteins.
thus, our finding of subtype-factors  would be the causes or indicators of disease  subtypes on genome level, which will be further studied in future.
CONCLUSIONS
pathway enrichment analysis is a useful tool in the study of biology or biomedicine, due to its functional screening on the well-known biological processes rather than single molecules.
the measurement of dysfunctions of pathways during a phenotype change, e.g., from normal to diseased, is the key issue when applying enrichment analysis for pathway or other functional gene set.
different from differentially expressed genes focused in previous methods, the genes with great differential expression variance are also attractive and important, which indicate another specific characteristic of a biological system in the change of phenotypes.
in the context of differential expression and expression variance, iea is proposed to identify the pathways full of degs and devgs simultaneously, rather than conventional approaches focusing on only degs.
the biological meaning of iea has obtained strong evidences by an evaluation scheme based on method comparison.
on the real datasets of disease samples, iea indeed specifically identify pathways containing degs and devgs, which are usually under-scored by other methods.
the map of pathways was further reconstructed based on the selected pathway crosstalks, and the module organization among pathways was also detected.
the topological structure of such network of pathways reveals the signalling pathways as upstream functions would be causes of disease, and the disease-relevant pathways as downstream functions would link to those upstream pathways by crosstalk.
in addition, some disease-relevant pathways or subtype-relevant pathways are well associated with clinical indices according to their devgs’ relative expression level, which are usually not observed from the raw expression profiles of pathway genes.
although many identified subtype-factors haven’t clinical evidences due to the limit in the clinical application, the iea actually show its ability to identify the risk of subtypes of genotype-phenotype associations.
those subtype-factors could help us in accurately realizing personal prevention or personal treatment  <cit> .
besides, the additional analysis results on colorectal cancer also support these conclusions .
totally, iea supplies a new way to carry on enrichment analysis in the context of differential expression and expression variance, and can easily expand to handle with the analysis in other more complicated context .
it is also necessary to expand iea to functional class scoring or pathway topology based approaches in future work.
additional files
additional file 1: table s <dig>  the 3-order approach-specific datasets corresponding to different methods based on prioritization performance.
additional file 2: table s <dig>  the prioritization performance about method comparison on approach-specific datasets .
additional file 3: table s <dig>  the prioritization performance about method comparison on approach-specific datasets .
additional file 4: table s <dig>  the sensitivity performance about method comparison on approach-specific datasets .
additional file 5: table s <dig>  the 3-order approach-specific datasets corresponding to different methods based on sensitivity performance.
additional file 6: table s <dig>   overlaps between pathway genes and disease associated genes on diabetes datasets.
additional file 7: table s <dig>  pathway ranking by different methods on diabetes datasets.
additional file 8: table s <dig>  pathway-pair as a crosstalk identified on diabetes datasets.
additional file 9: table s <dig>  association of each pair between a pathway and a clinic index  on diabetes datasets.
additional file 10: table s <dig>  association of each pair between a pathway and a clinic index  on diabetes datasets.
additional file 11: table s <dig>  association of each pair between a devg and a clinic index  on diabetes datasets.
additional file 12: table s <dig>  association of each pair between a devg and a clinic index  on diabetes datasets.
additional file 13:
si document - a case study on colorectal cancer and table a1-a <dig>
additional file 14: table a <dig>  pathways ranking for colorectal cancer.
additional file 15: table a <dig>  pathway-pair as a crosstalk .
abbreviations
degdifferentially expressed gene
devggene with great differential expression variance
ieaintegrative enrichment analysis
gseagene set enrichment analysis
rwrrandom walk with restart
slcsupervised-like clustering approach
pinprotein interaction network
competing interest
the authors declare that they have no competing interests.
authors’ contributions
tz and gl conceived of the study.
xy carried out the experiment and analysis.
xy and tz drafted the manuscript.
all authors read and approved the final manuscript.

 being a computational approach based on the prior knowledge, pathway enrichment analysis is widely used in the study of genotype-phenotype associations. biological pathway as a set of interactive genes produces particular cellular response/outcome by executing a series of functional cascades. it is curated by experts from wide range of science fields so that can supply more creditable functional details than general go module or network module. different from exploring the unknown or indeterminate functions by network module, pathway-centered analysis always makes an effort to capture the permutation of established functions in the change of phenotypes. as a key approach of pathway-centered analysis, the pathway enrichment analysis or well-known gene set enrichment analysis can identify dysregulated pathway by qualitatively measuring the changed status of a pathway. in the pathway enrichment analysis, the dysregulation of a pathway is the most important issue, and should be mathematically defined and measured well. it can estimate the conditional enrichment or status of a pathway, which is assumed to be associated with particular phenotypes. current researches generally use genes with significantly differential expressions or differential correlations to evaluate the extent of the dysregulation of a pathway. one kind of conventional method is evaluating the dysfunction of pathways in different conditions, such as fidepa, spia and ipeap. the other kind is using pathways to characterize individual samples, like corgs and pathifier. generally, all these methods focus on the genes with differential expression and their enrichments in pathways, which assume the samples are of good purity in genotype-phenotype association study. however, in the study of complicated phenotypes, e.g., cancer study, a relevant problem is the samples with the same disease phenotype might be full of different unknown subtypes due to disease heterogeneity. it is necessary to detect genes with new features observable in the complicated disease samples, and enhance the pathway enrichment analysis to be applicable in such previously unexpected situation. actually, there are new expression features extracted in recent studies, e.g., genes with differential expression variances. in the context of differential expression variance, it is still a challenge to measure the enrichment or status of a pathway. a solution to this problem can promote the efficiency of pathway enrichment analysis on genotype-phenotype association because it will consider more complete information about the expression changes of pathway genes. it can also provide new insights on the biological pathways by integrating additional expression and network features. in this work, we propose a multiple-label based enrichment analysis to detect such dysregulated pathways, which simultaneously takes into account the genes with differential expression and genes with differential expression variance together.fig. major differences between the measurements of dysregulated pathways used in conventional enrichment analysis and integrative enrichment analysis obviously, the hypothesis underlying iea is that the dysregulated pathways involved in disease heterogeneity would be full of degs and/or devgs. that means the identified pathways by iea would be disease pathways or their up-streams/down-streams. however, current methods in pathway enrichment analysis only expect to give high-rank to disease pathways. when iea identifies up-streams/down-streams of disease pathways, it further assistantly supplies a network of pathways to recover a global functional map and infer the associations among disease pathways and subtype-relevant pathways. noted, the biological meaning of the edge in such network of pathways is the pathway crosstalk, which is just an important biological mechanism or functional relationship among pathways. conventional researches tend to simply determine a pathway crosstalk by the overlapped genes in two pathways, which disregard the statistical significance of the genes and interactions involved in the pathway crosstalk. by contrast, degs and devgs in one pathway can be used as seeds, and further detected their interactive genes in the candidate crosstalking pathways by a random walk restart algorithm. the significance of a pathway crosstalk can be finally evaluated by the genes involved in this crosstalk as their enrichments in two pathways. based on the above concepts and mathematical models, a new pathway-centered analysis framework, the integrative enrichment analysis, is implemented as pathway enrichment score calculated by the hypergeometric test on differential genes; pathway crosstalk ranked by the random walk and hypergeometric test on rewired molecule networks; pathway-phenotype association and subtype-factors determined by devgs in pathways. according to a previously established evaluation scheme, we found that, in particular cohorts, a few target disease pathways can be significantly high-ranked by iea, which supplied the evidences of the deviation-based disease characteristics, and iea is more effective than other state-of-the-art methods in this condition. furthermore, by a proof-of-concept study, we shows the details of iea on analyzing real transcriptional data related to complex diseases, e.g., diabetes and colorectal cancer. iea indeed captures the previously under-estimated pathways full of devgs and degs. these newly identified dysregulated pathways would be heterogeneity-relevant pathways and are found to be significantly linked to disease pathways by estimated crosstalks. many candidate subtype-factors are also recognized as devgs or pathways associated with the risk of subtypes of genotype-phenotype associations. totally, iea supplies a new way of over-representation approach to carry on enrichment analysis in the complicate context of clinical application, and could be easily expanded to functional class scoring or pathway topology based approaches, which will be a necessary complementary and cooperative approach to conventional ones. the matlab scripts of the software named ieapackage and some alternative r scripts have been deposited in github and accessed in https://github.com/bluesky2009/integrative-enrichment-analysis. this software has been developed and tested in windows or windows and matlab or matlab methods generally, enrichment analysis includes three categories of methods: over-representation approach, functional class scoring and pathway topology based approaches. although these methods are all focusing on evaluating the phenotype-associated pathway, they would be based on different hypothesis. this work and the proof-of-concept study are based on the over-representation approach, which measures the dysregulation extent of a pathway according to the number of dysregulated genes in this pathway. traditional methods only evaluated the degs in a pathway; by contrast, iea evaluates the degs and devgs in a pathway. thus, the meaning of the statistic for the integration of iea is as completely as possible to measure the dysregulation extent of a pathway according to the number of dysregulated genes in this pathway, which have been well defined and introduced in follows. differential gene expression and differential expression variance given a gene x has expression profiles in control and case samples as x and x’ respectively, the expression variance of this gene in control and case condition are e2) and e2) respectively. here, u and u’ are average expressions of gene x in control and case samples respectively. then, the conventional criterion and measurement of genes with differential expression are: h0:ex=ex’;h0rejected; where x or x’ are the original/raw expression levels. noted, the differential expression includes up-regulation and down-regulation. except for these degs, the genes with differential expression variance are also discriminative features. the expression variance concerned features, e.g., bimodal gene expression, is already known as an important expression pattern in the control of a transition of biological systems, such as: disease development, cellular differentiation, and phase transition. however, the differential expression variance of genes has not been studied in a systematic way to the best of our knowledge, especially for its usage in the pathway enrichment analysis. the differential expression of genes, used in conventional enrichment analysis, requires the gene’s expressions under different conditions to distribute around different mean expression levels. by contrast, differential expression variance of genes can be defined as the genes’ deviations being significantly different under dissimilar conditions, such as: h0:ex‐u=ex’‐u’;h0rejected;andh0:ex=ex’;h0notrejected where x or x’ is the original expression level, |x-u| or |x’-u’| is the relative expression level. noted, the differential expression variance includes tight-regulation and relax-regulation. and importantly, as defined above, the devgs have excluded degs, or there is no overlap between devgs and degs in this work. that means, when one gene has both differential expression and differential expression variance, this gene is thought as deg in priority in order to be consistent with conventional analysis; and, of course, this kind of genes are worthy of deep research in future work. actually, given x or x’ satisfy normal distribution, |x-u| or |x’-u’| will be folded normal distribution, then the wilcoxon rank sum test instead of student’s t-test is used in the significance test of devgs. integrative enrichment analysis in the context of differential expression variance obviously, the conventional enrichment analysis limits to estimate the extent of differential expression rather than differential expression variance. when considering the contribution of devgs on pathway’s dysregulation, it is necessary to refine the conventional approach to take into account the degs and devgs together. naturally, an easiest strategy is to put degs and devgs together as the same dysregulated genes and use conventional hypergeometric test to obtain the p-value. however, this will disregard the respective distribution of degs and devgs in a target pathway and in the whole transcriptome. thus, we extended the hypergeometric test on two kinds of enriched genes simultaneously as bellows. our approach, noted as ht, still depends on the hypergeometric distribution and uses p-value to measure the dysregulation of a pathway in the context of differential expression variance. briefly seen in table  given there are expression data on total n genes, and x degs and x devgs selected respectively. for some pathway, k and k genes from pathway members have differential expression and differential expression variance respectively. then the significance of deregulated genes as degs or devgs enriched in this pathway can be estimated by formula this p-value also ranges from zero to one. the less the p-value is, the larger dysregulation extent the pathway has, when the significantly larger number of genes in this pathway show differential expression or differential expression variance. px1=k x2=k2=x1k1x2k2n−x1−x2y−k1−k2nypx1>k x2>k2=1−∑<i i2>∈0x1×0x2−k x1× and disease-relevant pathways. given several genes in a pathway as seeds, iea uses random walk to find their partner genes in the other pathway. in fact, random walk with restart is a well-known ranking algorithm for candidate gene prioritization. it supplies the probability of searching the random walker at nodes in the steady state, so that, it can give a measure of proximity between source nodes and other nodes in molecule network. let n be the adjacency matrix of a gene network with node set v and edge set e, in which the element nij equals one if e ∈ e represents the interaction between genes/nodes i and j), or zero otherwise. based on the topological structure of the gene network, the transition matrix t can be calculated. each element in the transition matrix is denoted as tij and represents the probability of transition from node i to node j. the value of tij can be given by one of two ways as follows, the first one is topology-weighted and the second one is correlation-weighted. tij=nijdi,ifeij∈e otherwise,wheredi=∑j∈vnij tij=wijnijwi,ifeij∈e otherwise,wherewi=∑j∈vwijnij the rwr algorithm updates the probability vectors by pk+1=1−λtpk+λp k> where t is the transition matrix and p is the initial probability vector with the sum of the probabilities as one. in p all the source nodes are assigned equal probabilities and other nodes are given zero. p∞ is obtained when the algorithm is convergent. if p∞ > p∞, node i is thought to be more proximate to source nodes than node j does. thus, a two-way rwr approach is proposed to search the genes involved in two interactive pathways and estimate their enrichment for evaluating the pathway crosstalk. the steps of two-way rwr include: for each pathway u, its degs and devgs are used as source nodes/genes, and rwr is used to rank the genes in known molecule network, e.g., protein association network collected from string database. in the high-ranked genes from above rwr analysis, the genes belonging to pathway v are the partner genes interactive with source genes. based on the sources genes and their partner genes, the enrichment of those interactive genes in pathways u and v can be evaluated by our ht approach, i.e., p-value in formula for every pathway, the analysis in steps and is repeated. then, given a pathway pair, it is a pathway crosstalk only when euv and evu are both significant. finally, the map of pathways consist of those selected pathway crosstalks, where a node represents a pathway and an edge represents a pathway crosstalk. screening subtype-factor of genotype-phenotype associations based on devgs and dysregulated pathways supplied by iea the second assistant down-stream analysis method of iea is to screen subtype-factors according to the available clinical indices. as stated above, iea focus on the devgs and their involved pathways, and these genes and pathways are thought as signatures of potential subtypes of heterogeneous samples. however, these hidden subtypes might have not been identified or formalized in clinics. to evaluate such new signatures or subtypes, one direct strategy is to measure the correlation between genetic signatures and clinical indices. if one signature is significantly related to some clinical index, the subtype represented by such signature would be medical meaningful as to be observable in clinics and this signature is also called as subtype-factor related to particular clinical index. the approach to identify such subtype-factors is described in bellows. for each pathway, its devgs are used to group case samples into two clusters, when the case samples have high varying expression compared to control samples. that means these genes have over-expression in one group of samples and under-expression in the other group of samples. this pathway would be a candidate subtype-factor when these two sample clusters are discriminative on some clinical index. on this condition, a clinical subtype of samples is thought to be related to a given clinical index, which is represented by a subtype-factor. obviously, the clinical subtype of a particular sample might be contributed by many subtype-factors. given a known phenotype, a few subtype-factors correlated with this phenotype can be found, although which just reveals only the tip of the iceberg for the subtypes of genotype-phenotype associations. particularly, different from conventional un-supervised clustering for subtype identification, a supervised-like clustering approach is proposed to identify subtype-factors on the level of pathways. firstly, the case samples can be grouped into two clusters according to their features’ values compared to those values of control samples: on each feature, one group of samples have larger values than controls meanwhile the other group of samples have less values than the same controls, or vice versa. that means, a hyperplane determined by a few control samples could separate the samples space into two sub-spaces, and case samples in each of two sub-spaces are grouped into one cluster. secondly, some clinical information of samples can be used to evaluate the potential subtype represented by such two clusters of case samples. if the clinical values of these two groups of samples have significant difference, a clinical subtype of genotype-phenotype association is identified and the corresponding pathway is a subtype-factor corresponding to the given clinical index. practically, the slc algorithm on a pathway is implemented as bellows: discrete the expressions of devgs of case samples into binary vector based on the values of controls: for a devg, if its expression value is larger than the mean of controls, it is one in the binary vector; otherwise, it is zero. clustering case samples based on the binary vectors by conventional methods as hierarchical clustering or k-means, which obtains two sample clusters. calculating the significance of difference between clinical indices among above two sample clusters. if the difference is significant, this pathway is identified as a subtype-factor of the association between the given pathway and clinical index.  the evaluation of biological meaning of iea by method comparison iea is proposed to evaluate dysregulated pathways by differential gene expression and differential expression variance together. differential expression variance has been reported as a new and important expression change during a phenotype change, e.g., diseases. in this work, the biological hypotheses underlying iea is that, the dysregulated pathways full of genes with differential expression variance would be subtype-relevant pathways. although subtype-relevant pathways for particular complex disease are unclear in current pathway databases, e.g., kegg, it is still able to investigate if prior-known disease pathways in kegg would be subtype-relevant and if iea can identify them. in the previous study of gene-set analysis, a comparison scheme has been built to evaluate the performances of different enrichment analysis methods based on multiple expression datasets about complex diseases. different from previous general comparison, we focus on the comparisons by approach-specific datasets, in order to mainly evaluate the biological meaning of iea. according to the comparison protocol, we ran total eight representative enrichment analysis methods on geo datasets with target pathways in kegg, and obtained the rank of target pathway estimated by each method on each dataset; then, for each dataset, we rank the eight methods according to their prioritization performance or sensitivity performance, and this dataset is assigned as a specific-data for the top-k methods; thus, all specific-data for one method can consist of k-order approach-specific dataset. generally, on one method’s approach-specific dataset, this method should have best or comparable performances than other methods, so that, the biological characteristics assumed by this given method would significantly displayed on these datasets. therefore, we can use this strategy to investigate the biological meaning of iea in real datasets. in bellows, we firstly summarize the biological hypothesis hold by different state-of-the-art enrichment analysis methods and their respective quantitative measurements, and then discuss the comparison between iea and others. plage: it assumes the activity of pathway rather than the expression of pathway genes determines the activated or inhibited status of pathways under different conditions; and the pathway activity is measured by an activity score as the weights of a metagene extracted from all pathway genes by svd. gsva: it proposes the change of pathway activity between control and case should be evaluated at the level of samples, e.g., considering the variation of pathway activity over a sample population; and the pathway activity is measured by so-called gsva score as a function of the expressions of genes inside and outside the pathway, and these scores are assessed similarly as gsea by using the kolmogorov-smirnov like random walk statistic. padog: it assumes that, if the genes highly specific to a given pathway occur differential expressions, the respective pathway would be truly relevant in that condition; thus, a new gene set score is calculated as the mean of absolute values of weighted moderated gene t-scores where the gene weights are designed to be large for the genes appearing in few pathways and small for genes that appear in many pathways. globaltest: it holds an assumption that, if a group of genes can be used to predict the clinical outcome, the expression patterns of such gene group must differ for dissimilar clinical outcomes; thus, it uses generalized linear model to give one p-value for a group of genes, not a p-value for each gene, which can be applied to estimate the enrichment of a given pathway. mrgse: it proposes that the high ranks of expression changes of genes can indicate the differential expression of a set of genes; and the enrichment score or the test statistic of a pathway is the mean rank of this gene set, i.e., the average of the ranks of t-statistics of pathway genes. gsa: it is similar to gsea, and proposes two improvements as the maximal average statistic for summarizing gene-sets, and restandardization for accurate enrichment inferences. ora: it takes into account the number of differentially expressed genes observed in a pathway as indicators of pathway states; generally, it uses a basic contingency table to test the association between the differential expression status of a gene and its membership in a given gene set, which can be measured by the p-value of a hypergeometric test. iea: it is proposed in this work to generally consider the contribution of expression variance in a dysregulated pathway; as one implementation, this work takes into account the number of degs and devgs observed in a pathway as indicators of pathway states; it is designed to test the association between the differential expression/differential expression variance status of a gene and their memberships in a given gene set, which can be measured by the p-value from proposed ht approach in this work. first of all, we can cluster the above eight approaches by their performances on all datasets to investigate the general association among different methods. as shown in figs.  and the similarity among any two methods is measured by four kinds of criterion: the first one is whether the ranks given by two methods on the same dataset are also the same; the second one is whether the ranks given by two methods have the same change tendency among different datasets; the third one is whether the p-values given by two methods on the same dataset are also the same; and the last one is whether the p-values given by two methods have the same change tendency among different datasets. obviously, gsa and padog are both based on conventional gsea, so that they are similar; the proposed iea is based on ora, thus, they also have similar performances on different datasets; plage and globaltest are closely clustered together, one reason is that they both estimate a score from all pathway genes rather than individual genes; in addition, mrgse and gsva are much different, and also different form other methods, which is possibly because they have specific design principles on the measurement of pathway dysfunctions, i.e., mrgse combines the t-statistics of individual pathway genes meanwhile gsva uses a score as a function of the expressions of genes inside and outside a pathway.fig. category of representative gene set analysis approaches based on clustering of prioritization performance. a method clustering based on euclidean distance of ranks of all pathways. b method clustering based on correlation distance of ranks of all pathways fig. category of representative gene set analysis approaches based on clustering of sensitivity performance. a method clustering based on euclidean distance of p-values of all pathways. b method clustering based on correlation distance of p-values of all pathways then, we directly grouped the datasets according to the performance of a given method, e.g., some datasets are included as k-order iea-specific datasets, only when the rank of iea performance compared to all methods are in the top-k on these datasets, where k is set in this study. to quantify the performance, sensitivity and prioritization are adopted as previously. in previous evaluation on these datasets, padog displays consistently comparable performance with other methods, meanwhile, plage, globaltest and mrgse have the best performances on some categorise of datasets, which already suggest the existence of approach preferences. thus, the categories of datasets induced in this work can indicate potential preferences of different methods on particular datasets. of course, on the preference-specific datasets, the preferred method should have the best performance; meanwhile, some other methods would have comparable performances. therefore, this new evaluation scheme can supply evidences for two hypotheses: one is that the expression features of a pathway measured by a given method actually have biological meaning, and are existing or observable in real datasets corresponding to particular phenotypes; the other one is that a method has significant preferences on a group of datasets, which is comparable or not by other methods. as seen in table  for prioritization performance, every method shows its preferences on particular datasets; in addition, the method comparison on more strict specificity as 2-order approach-specific datasets and more weak specificity as 4-order approach-specific datasets have also been done and reported in additional file 2: table s and additional file 3: table s combined these results together, it is easily to see that: iea, ora and padog have significant preferences, due to their highest performances than all other methods on their correspondingly preferred datasets; globaltest, gsva, plage and mrgse also have significant preferences on datasets, although sometimes padog could have comparable performances with them; besides, gsa is challenged by padog, because on the gsa-specific datasets, padog always achieves similar or better performances, but not vice versa.table the prioritization performance about method comparison on approach-specific datasets gsa padog iea mrgse ora globaltest gsva plage the performance of an approach on its specific dataset is highlighted in bold. and the performance of comparable approaches on some specific dataset is highlighted in bolditalic and seen in table  and additional file 4: table s for sensitivity performance, again, every method shows its preference on particular datasets. it seems that globaltest and plage have generally comparable performances with other methods according to their performances on many preferred datasets of other methods. even though, iea shows the best performance on the iea-specific datasets, and is better than ora on the ora-specific datasets. this fact strongly suggests that, iea actually can detect dysregulated pathways, and displays competitive performance than conventional methods in current evaluation scheme when those dysregulated pathways are just target pathways; besides, iea is realized based on the conventional ora, and improves ora on the sensitivity performance, which would just be contributed by considering the new feature genes as devgs in dysregulated pathways.table the sensitivity performance about method comparison on approach-specific datasets gsa padog iea mrgse ora globaltest gsva plage the performance of an approach on its specific dataset is highlighted in bold. and the performance of comparable approaches on some specific dataset is highlighted in bolditalic finally, above dataset-driven method comparison supplies new insights on the performance specificities of many gene-set approaches, especially for iea; and also supports the importance and biological meaning of dysregulated pathways identified by iea. we can draw following conclusions: although a few methods have consistent performances on many datasets, different method still have their preferences on the expression characteristics of dysregulated pathways, so that each method can achieve significantly better performances on their specific datasets rather than all datasets. especially, on the iea-specific datasets, iea indeed are the best one than other methods; and even on the ora-specific datasets, iea is better than ora on sensitivity performance. on the iea-specific datasets, the target pathways or disease pathways are possibly just the subtype-relevant pathways, so that, iea have competitive performance in the comparison scheme. therefore, the expression variance focused by iea is actually full of biological meaning, and will help iea to detect new dysregulated pathways, e.g., subtype-relevant pathways. in addition, the complex diseases concerned in iea-specific datasets actually already have reports about the existence of subtypes on genetic level, such as dilated cardiomyopathy, renal cancer, prostate cancer, colorectal cancer, and thyroid cancer. every method, or every method category, can actually capture particular dysregulated pathways. when the target pathway of a dataset just displays the expression characteristics focused by an approach, such approach would have better performance on this dataset. on some specific datasets preferred by other methods, iea should face two conditions: one is under the condition that the target pathway is not subtype-relevant pathway, and iea will be underestimated but has supplied a useful down-stream analysis to assistantly link the target pathway and potential subtype-relevant pathway identified; the other one is under the condition that the target pathway is a subtype-relevant pathway, and iea should be further enhanced by integrating expression variance with other pathway pattern, which is worthy of study in future. a proof-of-concept study of iea on transcriptional analysis of complex diseases iea has been applied to detect the biological malfunction of complex diseases on the pathway level rather than gene level. iea, as a pathway-centred analysis approach, not only supplies the conventional pathway enrichment analysis but also extracts divergent pathway associations, such as: pathway & disease genes, pathway & pathway, and subtypes of pathway & clinic. firstly, it is the data pre-procession. data needed in iea have been prepared from the public resources: the gene expression data of human islets from non-diabetic and diabetic were downloaded through geo; there are two datasets, the main dataset gse contains samples from non-diabetic and diabetic with genes, and the replicate dataset gse contains samples from non-diabetic and diabetic with genes; the gene lists of kegg pathways are obtained from gsea package; the human protein interaction network are extracted from string database with confidence score no less than; four clinical indices are also obtained from the supplementary of original study as sex, age, bmi, and hba1c; meanwhile, diabetes associated genes are searched from genecards database. secondly, it is the main step of iea. different scores of pathway enrichment are calculated: the conventional score as p-values of hypergeometric distribution of degs in a pathway; the conventional score as p-values of estimated pathway enrichment; the new score as p-values of hypergeometric distribution of differential genes in a pathway calculated by the proposed ht approach. for ora or iea, the thresholds of p-value of significance test on degs or devgs are both set as, and adopted a pervious strategy to select those feature genes: 1) select all genes with fdr adjusted p-values no more than; 2) if the genes selected are less than re-select all genes with p-values no more than and fold-change no less than; 3) if the genes selected are still less than directly use the top 1 % of genes ranked by p-values from least to largest. thirdly, it is one assistant step of iea. pathway crosstalks are evaluated by two-way rwr approach. the interactions selected from pin consist of differential network, where the selected interactions have significant correlation difference between diabetic and non-diabetic groups. on this differential network, two-way rwr approach is used to find the pathway crosstalks. the most significant pathway crosstalks consist of the map of pathways. besides, the enrichments of pathway genes or disease-associated genes in the high-ranked genes of rwr are also analyzed and evaluated by auc, which support the efficiency of rwr on pathway-related analysis. fourthly, it is the other assistant step of iea. the devgs in each pathway are used to group samples in two clusters by slc approach. for each clinical index from sex, age, bmi, and hba1c, its subtype-factors are identified. as comparisons, the all genes of each pathway are also used to directly group samples to evaluate the significance of detected subtype-factors. obviously, the above analysis routine can be applied on any other dataset of samples with different kinds of complex diseases. noted, our proof-of-concept study is to combine a group of genes with differential expression and a group of genes with differential expression variance. we don’t select the genes with high variances, but the genes whose expression variances can distinguish different conditions/phenotypes. if the expression variance of a gene can classify samples with different phenotypes well, the dominant component of this gene’s variance could be biological variance. obviously, the technique variance should have no such discrimination. we have also checked the correlation between the variance of each devg and each clinical index. many genes’ expression variances even have significant correlation with clinics, which would not be caused by technique variance too. our experiment on microarray is an application of iea to support the idea of combining deg and devg, and the results show our method’s efficiency. indeed, many approaches are still proposing to improve the selection of conventional deg or even new devg. the removal of technique variance will improve the selection of devg and final iea, which is our future work. diabetes associated genes on pathways first of all, we investigated the overlap between prior-known disease genes and pathway genes. many pathways are full of diabetes associated genes, which means pathways could have great changes during disease development and progression. these pathways would be causes or outcomes of the disease. although iea pays attention to the identification of dysfunctional pathways rather than discrimination of causal pathways, as introduced in follows, the map of pathways can further complementally supply some clues of the causal roles of pathways at the level of network of networks. then on the main dataset gse, we have obtained feature genes as summarized in table  there are degs and devgs selected by iea, many of them are also detected on the replicate dataset gse. the genes of degs are disease genes, and genes are pathway genes; meanwhile, the genes of devgs are disease genes and genes are pathway genes. obviously, there are many disease-informative or function-informative genes disregarded in conventional analysis, and iea can capture these genes and estimate their effects in the dysfunction of pathways.table the statistic on degs, devgs and their overlapping with pathway or disease genes adeg points genes with differential expression; devg points genes with differential expression variance; pg_deg points the pathway genes in degs, i.e., the overlaps between pathway genes and degs; dg_deg points the disease genes in degs, i.e., the overlaps between disease genes and degs; pg_devg points the pathway genes in devgs, i.e., the overlaps between pathway genes and devgs; dg_devg points the disease genes in devgs, i.e., the overlaps between disease genes and devgs; degup and degdown point genes with up-regulation and down-regulation respectively; devgup and devgdown point genes with relax-regulation and tight-regulation respectively furthermore, in degs, there are gene up-regulated in disease state and genes down-regulated. meanwhile, there are genes tight-regulated in disease condition and genes relax-regulated. the examples of such four expression patterns are shown in fig.  myc, known as a cancer oncogene, is also reported to be altered in diabetes. seeing fig. 4a, myc is indeed up-regulated in the diabetes state. insulin is known as a main cause of diabetes, and its an isoform as ins-igf actually has down-regulation when diabetes occurs. hoxd as a gene in the homeobox family encoding a highly conserved family of transcription factors, has an important role in the morphogenesis. it has tight-regulation in diabetes condition, so that, it may participate in the accurately regulation of biological processes associated to diabetes. by contrast, rexo known as transcription elongation factor b polypeptide 3-binding protein would be a cofactor involved in gene regulation. this gene shows relax-regulation in diabetes condition. thus, rexo would be a cause or indicator of some subtypes of diabetes. indeed, the original study has supplied four clinical indices, we found age is mostly related to the subtypes or sample clusters of diabetes determined by rexo.fig. the cases of differential expression patterns of deg and devg. for genes myc, ins-igf hoxd and rexo - give their expression profiles on dataset gse41762; - give their expression distribution on dataset gse41762; - give their expression profiles on dataset gse38642; - give their expression distribution on dataset gse38642 more importantly, devgs don’t mean no-changes. as the key point of our model and method, devg means a gene would denote activation of a signalling pathway in a group of samples, meanwhile, inactivation of this signalling pathway in another group of samples. this would be a main cause of heterogeneous samples. the biological mechanism underlying this phenomenon would be the switch of pathway activation. if based on other kinds of enrichment analysis framework, it is possible to discuss the activation, inactivation, or activation-switch of a signalling pathway, which will be studied in our future work. dysregulated pathways identified to capture degs and devgs simultaneously in conventional analysis, the genes with differential expression are focused; now, the genes with differential expression variance are also attractive. in the context of differential expression variance, the dysregulated pathways are expected to have as many degs & devgs as possible, which can be captured by over-representation approach like the proposed iea. for evaluation, a gene-distribution graph is further introduced to show the percentages of degs and devgs respectively for each pathway. in fig.  a pathway is represented by a point whose indices in axis are the percentages of degs and devgs in this pathway respectively. obviously, the pathways full of degs and devgs tend to locate at the right-up conner of such gene-distribution graph.fig. the gene-distribution graph of pathways identified by different methods. a-c the pathways high-ranked by ora, iea and gsea respectively, whose results on dataset gse d-f the pathways high-ranked by ora, iea and gsea respectively, whose results on dataset gse38642 to validate the pathways identified by iea and compare to other methods ora and gsea, we firstly calculated the percentages of degs and devgs of each kegg pathway and plot them on the gene-distribution graph. then we calculated the dysfunctional score and rank all kegg pathways, and found the top- selected pathway in graph. obviously, the pathways high-ranked by ora have least scores located at the right-bottom of gene-distribution graph. it means the pathways selected by ora are full of degs rather than devgs. by contrast, in the iea, the detected pathways are full of degs and devgs respectively and significantly. even more, the dysregulated pathways detected by well-known gsea on these datasets shown weak performance on the identification of pathways full of devgs. thus, iea indeed can effectively detect the pathways under-scored in conventional analysis, and these dysfunctional pathways would be disease-relevant or subtype-relevant. this conclusion is also supported by the similar results from the analysis on replicate dataset. noted, the pathways identified by different methods can be significantly observed in the analysis on the replicated dataset. the issue of the consensus of pathway identification is not discussed more here, although some other studies have worked to improve the robustness of dysregulated pathway identification by integrating other prior information. in the high-ranked pathways identified by iea, many pathways are actually full of devgs. these devgs could be further associated to potential subtypes of samples, which are discussed in follows. besides, some of these pathways under-scored by other methods indeed have been reported to be altered in the disease state. for examples,'kegg hematopoietic cell lineage'. diabetes is known to compromise the function of the bone marrow, and diabetic complications mainly including macrovascular events might be from the dysfunctional bm-derived hematopoietic cells.'kegg cytokine-cytokine receptor interaction'. cytokines regulate inflammatory and immune responses, which play important roles in the pathogenesis of diabetes and its microvascular complications. the functional variations of cytokines and their receptors can benefit the prediction of the susceptibility and progression to diabetic nephropathy. as the potential pathogenic mediators in dn, cytokines might provide new potential therapeutic agents for disease treatments. dysregulated pathway-crosstalk identified to reveal the interactive map and module among pathways different from conventional strategy to use the overlapping genes as pathway crosstalks, two-way rwr is assistantly applied to find the interactive genes between any two pathways. rwr is previously used in the rank of disease genes, which holds the assumption that the candidate pathogen genes are more proximate to the known disease genes than randomly selected genes. this assumption is expected to be hold for pathway genes too. we first evaluated the pathway genes possibly selected by rwr,which uses the identified degs & devgs in a pathway as seeds. in the high-scored genes by rwr, there is a significant amount of pathway genes, which is the same as disease genes ranked by rwr. in the evaluation in table  we have set two kinds of control experiments. one is the prior-known network used, i.e., two sources as string and hprd are both applied; the other one is the feature genes used in the given background network, i.e., three kinds of feature genes are respectively used to calculate the auc values to evaluate the efficiency of selecting/ranking pathway genes or disease-associated genes. depending on these experiments, we can find:table the auc of different rank lists for pathway genes afor feature genes like all genes and degs, they don’t use network information, so that, they have the same auc values on the same dataset although different network used afor feature genes like all genes and degs, they don’t use network information, so that, they have the same auc values on the same dataset although different network used as a control, when all genes are ranked according to p-values, its auc is low. and when the selected degs are ranked according to p-values, the auc is similar to that of all genes. meanwhile, when the genes from two-way rwr excluding seeds are ranked according to proximity values, the auc achieves highest, which support again rwr-based approach is effective to capture interactive phenotypic genes as pathway genes or disease-associated genes. thus, two-way rwr is effective to mimic the pathway crosstalk and construct the associations among pathways, which is obviously consistent in multiple control experiments by using different prior-known network, or different datasets, or even different ranked/selected feature genes. the protein association network rather than protein physical network would be more efficient to lead the two-way rwr to link seed genes to pathway genes or disease-associated genes. there are at least two reasons for this result: one is that the known protein physical network is still greatly incomplete, by contrast protein association network would supply additional predicted interactions with high confidence; the other one is that, except for direct interactions between pathways, protein association network would cover much more indirect interactions or long-term interactions, which would mimic the pathway crosstalk well. in all, protein association network is efficient to detect the associations among pathways, whose usage to accurately predict physical interaction of pathways is out of this work and would be a future study. thus, rwr is actually effective to detect the genes interactive within known pathway genes. then by two-way rwr, we can find the interactive genes from two pathways, and select any pathway-pair as a crosstalk significantly. all the crosstalks connect the known pathways as a map, where each crosstalk is an edge and a pathway is a node. this map of pathways is a network of networks, rather than the original background network of separate molecules. in the map of pathways, the modules of interactive pathways can be detected, where a module represents a group of closely inter-connected pathways. in these pathway modules, two modules are obviously related to diabetes. one module is the group of signalling pathways. signalling pathways are known as the up-stream functions in the cascades of signals, so that, they have great possibility to be the causes of the dysfunction of down-stream functions, e.g., diabetes pathways. the other module is just the group of pathways concerning glycolysis, sugar metabolism, glycosaminoglycan and diabetes, which seems to be extremely a core pathway module of diabetes. noted, the type ii diabetes pathway would have significant crosstalk with pentose phosphate pathway as shown in fig.  as reported, the pentose phosphate pathway is widely activated in diabetes and its complications, thus this pathway would be important to understand the risk of diabetes diagnosis and treatment in clinical application. obviously, the pentose phosphate pathway have no significance on the selection by iea or other methods, and actually, its importance is reflected from the topological structure of the map of pathways. these facts reveal: conventional approaches usually focus on single pathways, so that, they can sometimes capture the disease associated pathways relevant to particular phenotypes but can’t distinguish or underestimate the potential causal relationship among pathways; the proposed iea supplies the map of pathways to reflect the functional organization of pathways, and disclose the key modules of pathways, such as the upstream pathways related to signalling pathways and the downstream pathways associated to diseases; on the map of pathways, those pathways full of devgs and degs tend to interact with disease pathways, indicating the determinant of subtypes are actually also the potential determinants of diseases. totally, the map of pathways supplies us a new viewpoint of functional organization at the level of network of networks.fig. the main topological structure of the map of pathways reconstrcucted on diabetes datasets. module is full of signalling pathways, and module is full of disease/diabetes associated pathways dysregulated pathways associated to clinical indicator as subtype-factors the pathways full of degs and devgs are specially selected by iea, in which the devgs might be the cause of potential subtypes of samples. compared to the mean values of genes in case, the samples in control can be divided into two groups: one group of samples have larger gene expressions than that mean value; and the other group of samples have fewer values, or vice versa. thus, we have applied two strategies to cluster the samples in control or case group, which can associate one pathway to some clinical index. the first common strategy is using the expressions of whole genes in a pathway to cluster samples into two clusters, and test the significance of these two groups of samples on one clinical indicator. the second strategy is using the discrete value of devgs in the same pathway to group samples and measuring their relationships with clinical indices. additional file 9: table s gives the p value of the association of each pair between a pathway and a clinic index for normal samples, and additional file 10: table s gives those for diabetes. generally, slc tends to discover more significant potential subtypes of samples corresponding to particular clinic index. thus, devgs actually have more power to identify the subtypes of genotype-phenotype associations than conventional approaches based on differential expression only.fig. the summary of dysregulated pathways identified as subtype-factors corresponding to four clinical indices. for example, there are pathways identified by iea to be correlated with hba values, but, there are only pathways recognized by conventional strategy. similarly, for other clinical indices, iea can also identify many subtype-factors obviously, one pathway can associate to multiple clinical indices, and one clinical index can also relate to multiple pathways. ‘kegg wnt signaling pathway’ is consistently related to bmi in the analysis of duplicate datasets. in fact, this pathway already has a few evidences on their associations with body weight. ‘kegg cell cycle’ is found to possibly associate with age, sex, bmi or hba, which is known as a common factor in disease development and progression. although each devg can separately determine some potential subtype-factor, here, we only discuss the determination of devg combination at the level of pathway. in addition, the association of each pair between a gene and a clinical index are listed in additional file 11: table s and additional file 12: table s for normal and diseased samples respectively. noted, subtype is a potential biological explanations of devg, so that, we have evaluated the possible subtype-factors rather than subtypes by the association between devg and known clinical index. in the statistical analysis, we have analyzed age, sex, bmi, and even hba1c which evaluates the risk of diabetes. thus, this clinical information can help us to find subtype of genotype-phenotype associations as stated above. of course, in breast cancer, there are some well-known subtypes determined on genotypes. to the best of our knowledge, in diabetes, the subtypes as t1d, t2d, gestational diabetes, surgically induced diabetes, chemically induced diabetes, are not defined by one or two genes/proteins. thus, our finding of subtype-factors would be the causes or indicators of disease subtypes on genome level, which will be further studied in future.  pathway enrichment analysis is a useful tool in the study of biology or biomedicine, due to its functional screening on the well-known biological processes rather than single molecules. the measurement of dysfunctions of pathways during a phenotype change, e.g., from normal to diseased, is the key issue when applying enrichment analysis for pathway or other functional gene set. different from differentially expressed genes focused in previous methods, the genes with great differential expression variance are also attractive and important, which indicate another specific characteristic of a biological system in the change of phenotypes. in the context of differential expression and expression variance, iea is proposed to identify the pathways full of degs and devgs simultaneously, rather than conventional approaches focusing on only degs. the biological meaning of iea has obtained strong evidences by an evaluation scheme based on method comparison. on the real datasets of disease samples, iea indeed specifically identify pathways containing degs and devgs, which are usually under-scored by other methods. the map of pathways was further reconstructed based on the selected pathway crosstalks, and the module organization among pathways was also detected. the topological structure of such network of pathways reveals the signalling pathways as upstream functions would be causes of disease, and the disease-relevant pathways as downstream functions would link to those upstream pathways by crosstalk. in addition, some disease-relevant pathways or subtype-relevant pathways are well associated with clinical indices according to their devgs’ relative expression level, which are usually not observed from the raw expression profiles of pathway genes. although many identified subtype-factors haven’t clinical evidences due to the limit in the clinical application, the iea actually show its ability to identify the risk of subtypes of genotype-phenotype associations. those subtype-factors could help us in accurately realizing personal prevention or personal treatment. besides, the additional analysis results on colorectal cancer also support these conclusions. totally, iea supplies a new way to carry on enrichment analysis in the context of differential expression and expression variance, and can easily expand to handle with the analysis in other more complicated context. it is also necessary to expand iea to functional class scoring or pathway topology based approaches in future work. additional files additional file 1: table s the 3-order approach-specific datasets corresponding to different methods based on prioritization performance. additional file 2: table s the prioritization performance about method comparison on approach-specific datasets. additional file 3: table s the prioritization performance about method comparison on approach-specific datasets. additional file 4: table s the sensitivity performance about method comparison on approach-specific datasets. additional file 5: table s the 3-order approach-specific datasets corresponding to different methods based on sensitivity performance. additional file 6: table s overlaps between pathway genes and disease associated genes on diabetes datasets. additional file 7: table s pathway ranking by different methods on diabetes datasets. additional file 8: table s pathway-pair as a crosstalk identified on diabetes datasets. additional file 9: table s association of each pair between a pathway and a clinic index on diabetes datasets. additional file 10: table s association of each pair between a pathway and a clinic index on diabetes datasets. additional file 11: table s association of each pair between a devg and a clinic index on diabetes datasets. additional file 12: table s association of each pair between a devg and a clinic index on diabetes datasets. additional file 13: si document - a case study on colorectal cancer and table a1-a additional file 14: table a pathways ranking for colorectal cancer. additional file 15: table a pathway-pair as a crosstalk. abbreviations degdifferentially expressed gene devggene with great differential expression variance ieaintegrative enrichment analysis gseagene set enrichment analysis rwrrandom walk with restart slcsupervised-like clustering approach pinprotein interaction network competing interest the authors declare that they have no competing interests. authors’ contributions tz and gl conceived of the study. xy carried out the experiment and analysis. xy and tz drafted the manuscript. all authors read and approved the final manuscript. 
however, the disease samples are usually heterogeneous, so that, the genes with great differential expression variance are becoming attractive and important to indicate the specific state of a biological system. in the context of differential expression variance, it is still a challenge to measure the enrichment or status of a pathway. one is to infer the association among identified dysregulated pathways and expected target pathways by estimating pathway crosstalks. furthermore, we present a proof-of-concept study on diabetes to indicate: iea rather than conventional ora or gsea can capture the under-estimated dysregulated pathways full of devgs and degs; these newly identified pathways could be significantly linked to prior-known disease pathways by estimated crosstalks; and many candidate subtype-factors recognized by iea also have significant relation with the risk of subtypes of genotype-phenotype associations. totally, iea supplies a new tool to carry on enrichment analysis in the complicate context of clinical application, as a necessary complementary and cooperative approach to conventional ones. the differentially expressed genes are widely focused in conventional analysis, which is based on the great purity of samples. pathway enrichment analysis is a useful tool to study biology and biomedicine, due to its functional screening on well-defined biological procedures rather than separate molecules. based on a previously established evaluation scheme, we found that, in particular cohorts, a few target disease pathways can be significantly high-ranked by iea, which is more effective than other state-of-the-art methods. the measurement of malfunctions of pathways with a phenotype change, e.g., from normal to diseased, is the key issue when applying enrichment analysis on a pathway. to address this issue, we proposed integrative enrichment analysis based on a novel enrichment measurement.
 pathway enrichment analysis is a useful tool to study biology and biomedicine, due to its functional screening on well-defined biological procedures rather than separate molecules. the measurement of malfunctions of pathways with a phenotype change, e.g., from normal to diseased, is the key issue when applying enrichment analysis on a pathway. the differentially expressed genes are widely focused in conventional analysis, which is based on the great purity of samples. however, the disease samples are usually heterogeneous, so that, the genes with great differential expression variance are becoming attractive and important to indicate the specific state of a biological system. in the context of differential expression variance, it is still a challenge to measure the enrichment or status of a pathway. to address this issue, we proposed integrative enrichment analysis based on a novel enrichment measurement.  the main competitive ability of iea is to identify dysregulated pathways containing degs and devgs simultaneously, which are usually under-scored by other methods. next, iea provides two additional assistant approaches to investigate such dysregulated pathways. one is to infer the association among identified dysregulated pathways and expected target pathways by estimating pathway crosstalks. the other one is to recognize subtype-factors as dysregulated pathways associated to particular clinical indices according to the devgs’ relative expressions rather than conventional raw expressions. based on a previously established evaluation scheme, we found that, in particular cohorts, a few target disease pathways can be significantly high-ranked by iea, which is more effective than other state-of-the-art methods. furthermore, we present a proof-of-concept study on diabetes to indicate: iea rather than conventional ora or gsea can capture the under-estimated dysregulated pathways full of devgs and degs; these newly identified pathways could be significantly linked to prior-known disease pathways by estimated crosstalks; and many candidate subtype-factors recognized by iea also have significant relation with the risk of subtypes of genotype-phenotype associations.  totally, iea supplies a new tool to carry on enrichment analysis in the complicate context of clinical application, as a necessary complementary and cooperative approach to conventional ones. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. issue-copyright-statement© the author 2015 