BACKGROUND
dna methylation is a chemical modification of cytosine bases, which is critical for cellular differentiation, cell development and disease progression .
for example, dna methylation directly inhibits the binding of transcription factors  <cit> , and methylation aberrations either predispose to or result in disease progression  <cit> .
with biotechnological advancements, dna methylation is also considered a biomarker of epigenome analysis  <cit> .
multiple platforms, including reduced representation bisulfite sequencing  <cit> , whole genome shotgun bisulfite sequencing  <cit> , methylc-seq  <cit> , and target capture bisulfite sequencing  <cit>  and dna methylation beadarray, have been developed to generate genome-wide dna methylation data.
high-throughput technologies have generated large-scale genome-wide dna methylation profiles for various cancers and cell lines, providing great opportunities for revealing the epigenetic mechanisms.
various approaches have been proposed on the basis of methylation profiles to extract dna methylation patterns.
for example, fleischer et al.
<cit>  identified  <dig> cpg probes associated with the survival time of breast cancer patients.
hinoue et al.
<cit>  recognized four distinct subgroups in colorectal cancer by analyzing large-scale genome-wide dna methylation profiles.
yang et al.
<cit>  discovered the comethylation modules across  <dig> cell lines by using a weighted comethylation network.
varley et al.
<cit>  identified dynamic dna methylation patterns across  <dig> human cell lines.
gevaert et al.
<cit>  designed methylmix to obtain  <dig> pancancer clusters which reveal a novel similarity across various cancers.
further details have also been provided  <cit> .
although dna methylation patterns have been extensively investigate, open questions have yet to be answered.
for example, it is acknowledged that dna methylation aberrations cause diseases by mediating gene expression  <cit> .
however, how the dna methylation regulates gene expression remains unknown.
epigenetic modules is critical for revealing the epigenetic mechanisms of cancers.
thus, algorithms for functional epigenetic modules should be established by simultaneously analyzing methylation and gene expression data.
however, the design of algorithms for the integrative analysis of dna methylation and gene expression is highly nontrivial because of two reasons.
first, integrative analysis requires large-scale sample matched methylation and gene expression profiles.
second, dna methylation patterns are difficult to determine because the relationship between methylation and expression is unknown.
for instance, it is acknowledged that the correlation between promoter methylation and gene expression is negative.
however, the recent evidence indicates that the positive correlation is also presented  <cit> .
regarding the first issue, the consortium, such as tcga, has generated sample matched dna methylation and gene expression data for various cancers  <cit> .
thus, functional epigenetic module can be identified.
regarding the second issue, researchers developed many algorithms, such as the epimod algorithms  <cit> , based on the gene comethylation network by using methylation data and searching modules in the network.
however, these algorithms are limited because they are solely based on either methylation or gene expression.
consequently, these algorithms fail to obtain epigenetic modules precisely.
to address this problem, researchers developed the fem algorithm  <cit>  that integrates dna methylation and gene expression into a protein interaction network by assuming the anti-correlation between dna methylation and gene expression.
the algorithm successfully identifies a novel epigenetically deregulated hotspot and methylated gene modules.
this finding indicates the superiority of integrative analysis in terms of revealing methylated gene modules.
current algorithms are limited by heterogeneous data integrated into a network, and thus require the correlation between dna methylation and gene expression.
however, the correlation between heterogeneous data is unknown and difficult to determine.
recently, multiple networks have been widely used to characterize the complex biological patterns.
for example, ma et al.
used multiple networks to characterize the dynamic modules  <cit> .
cantini et al.
proposed a novel multi-network-based strategy to integrate different layers of genomic information and developed algorithm to identify cancer driver genes  <cit> .
didier et al.
assessed aggregation, consensus and multiplex-modularity approaches to detect communities from multiple networks  <cit> .
these algorithms indicate that characterizing biological patterns in multiple networks is more accurate than those in a single network.
therefore, we propose an alternative method that characterizes the epigenetic modules by using multiple networks.
we first construct the differential comethylation and coexpression networks, and define the common modules within multiple networks as the epigenetic modules .
the proposed strategy is advantageous because it avoids specifying the correlation between heterogeneous data and the accuracy in the identification of epigenetic modules is improved.
we also develop the emdn algorithm by simultaneously analyzing multiple networks .
experimental results demonstrate that the proposed algorithm is more accurate than other algorithms.
we further demonstrate that epigenetic modules enriched by the known pathways, serve as biomarkers to predict breast cancer subtypes and survival time of patients.
the proposed model and algorithm provide an effective way for the integrative analysis of dna methylation and gene expression.
fig.
<dig> flowchart of the proposed algorithm.
a the common modules in both gene expression and methylation networks are defined as the functional epigenetic modules,  emdn consists of two major components: differential network construction and functional epigenetic module discovery.
the differential network construction consists of two steps: first, the coexpression  network is constructed using the gene expression  profiles, then the p-values of gene expression  between tumor and non-tumor samples are incorporated into the coexpression  network to develop the differential networks.
finally, the algorithm discovers the common modules in multiple differential networks
methods
the emdn consists of two components: network construction and epigenetic module discovery .
step 1: differential network construction
the development of a gene differential coexpression network involves two steps: 1) the development of a binary coexpression network and 2) edge weight assignment based on differential gene expression between normal and tumor conditions.
we assume that the gene expression data represents normalized estimates of gene expression intensity and are summarized at the gene level, which includes rna-seq or expression data generated by using illumina beadchip or affymatrix arrays.
the binary gene coexpression network is developed where the edge weight on a pair of genes is the absolute pearson coefficient of the gene expression profiles.
to remove the indirect correlation due to a third gene, w the first order partial pearson correlation coefficient is used  <cit> .
finally, edges whose weights are equal or greater than the pre-defined threshold δ are selected.
using the limma package  <cit> , the p-values of the gene expression difference between normal and tumor samples for each gene are obtained.
the p-values are corrected by benjamini-hochberg   <cit> .
then, the gene differential coexpression network is developed by using the p-values.
specifically, the weight w
ij on edge  in the differential network is defined as
<dig> wij=logpi+logpj1/22maxllogpl1/2if|cor|≥δ,0if|cor|<δ,
where p
i is the p-value of the differential expression for gene i, and c
o
r is the value of pearson correlation between the i-th and j-th gene based on the expression profiles.
the genes that are co-expressed and significantly differentially expressed are assumed to have assigned heavy weights.
in this study, the δ is  <dig> .
this value leads to the maximal number of genes connected in all the networks.
analogously, the differential comethylation network is developed using the gene methylation profiles.
step 2: discovering the epigenetic modules
to infer the epigenetic modules, emdn is based our recently developed m-module algorithm that is designed for identifying common modules in multiple molecular interaction networks  <cit> .
emdn consists of three steps: seed prioritization, module search by seed expansion and refinement of candidate modules.
the seed prioritization ranks genes in multiple networks by using the topological feature of the gene in the networks.
specifically, for each network g
k with the weighted adjacency matrix w
k=n×n, a function f
k:v→r is developed, where f
k denotes the importance of the i-th gene in the corresponding network.
the function is defined as
<dig> fk=αak′fk+y,
where ak′fk denotes the topological importance of nodes and y is the vector for the prior information.
the parameter α controls the relative weight of the topological importance and prior knowledge.
ak′ is normalized adjacency matrix, i.e., ak′=dk−1/2wkdk−1/ <dig> with d
k=d
i
a
g. we use the following iteration-based algorithm to obtain f
k:
<dig> fk=αak′fk+y,
where t denotes the iteration, and fk <cit> = <dig>  no prior information is used.
thus, we set y= <dig>  the iteration is terminated if there is no change between f
and f
.
usually, only  <dig> iterations are observed.
we calculate the z-score for f
k. after ranking the genes in all networks, the ranks of genes are obtained in multiple networks, denoted by r=[f
<dig> f
<dig> …,f
m].
the final gene ranking is obtained by using the z-score of genes in multiple networks, i.e., the row sum of r.
starting with each seed, the module search step iteratively adds genes whose addition causes the maximum decrease in the graph entropy-based objective function until no decrease in the objective function is observed.
given a module c, the entropy function for network g
k is defined as
<dig> hk=−pklogpk−log,
where pk=∑i,j∈cwijk/∑i∈cwijk.
the entropy function for c across all networks is defined as
<dig> h=∑k=12hk/|c|.
the module search procedure expands the module by using eq.
.
each time the gene whose addition leads to the maximum decrease of h is included in c until no gene can improve h.
during the refinement step, candidate modules whose sizes are smaller than  <dig> are removed.
to merge overlapping function modules, two modules whose jaccard index is greater than  <dig>  are merged.
data
we download gene expression and dna methylation data for breast cancer from tcga database.
specifically,  <dig> samples with matched level- <dig> illumina 450k methylation data and hiseq rsem gene-normalized rna-seq data are obtained with  <dig> tumor and  <dig> normal samples.
for methylation data, the β signal of the probe is used, which is calculated as the methylated signal divided by the sum of the methylated and unmethylated signal.
for the rna-seq data, the reads per kilobase of exon model per million mapped reads  is used.
the clinical information is also obtained from tcga.
in all these datasets, probes with more than 30% missing values are removed, and probes with less than 30% missing values are imputed using the r package pamr  <cit> .
gene methylation profiles
to assign dna methylation to a given gene , we follow the strategy in the ref.
<cit> .
specifically, the average value of the probes mapping within  <dig> bp of the transcription start site  is assigned to the gene.
if no probes mapped within  <dig> bp of the tss, we use the average value of probes mapping to the 1st exon of the gene.
if such probes are also not available, we use the average value of probes mapping within  <dig> bp of the tss.
eigengenes of modules
the eigengene of a module is defined as the first principal component based on singular value decomposition   <cit> .
in details, the gene expression matrix of a given module is denoted by x= where the index i corresponds to the module genes and the index j corresponds to the samples.
the singular value decomposition of x is denoted by
<dig> x=udvt,
where the columns of the matrices u and v are the left- and right-singular vectors, respectively.
the first column of v is the module expression eigengene.
similarly, we obtain the module methylation eigengene by the gene methylation profiles.
correlation of module methylation and expression
the correlation between module methylation and expression is defined as the pearson correlation coefficient between module methylation and methylation eigengenes.
the correlation is significant when bh correction is used with the cutoff  <dig> .
if the sign of pearson coefficient is greater than  <dig>  a positive correlation is observed, and negative if otherwise.
survival analysis
we use the function coxph   <cit>  to implement the cox proportional hazard model to analyze the association of methylation profile of each epigenetic module with the patient survival.
we use the prognostic index to generate high- and low-risk patient groups
<dig> indexi=∑c=1kβcxci
where k is the number of cancer-specific modules, β
c is the regression coefficient of the cox proportional hazard model for the c-th module and x
ci is the average methylation level of genes within the c-th module in the i-th patient.
patients are grouped into high- and low-risk groups based on the median of prognostic index.
the survival difference between these two groups of patients is obtained by using the kaplan-meier estimator and log-rank method.
hand <dig> module and simulation
as the true module, we select the hand <dig> module , since the biological and clinical significance of the driver gene, hand <dig>  has been extensively validated  <cit> .
to fully assess the sensitivity and specificity of the proposed method, the simulation model of hand <dig> module is also adapted  <cit> , in which it simulates statistics of differential methylation and differential expression on the protein interaction network.
the model bootstrapped statistics for the member genes of this module to come from the top and lower 5% statistics quartiles, with the statistics of the rest of the network genes bootstrapped from the middle 90% portion.
for each simulation run, the accuracy is recorded.
statistical significance of modules
the statistical significance of modules is based on the null score distribution of the random modules generated using randomized networks.
each network is completely randomized  <dig> times by degree-preserved edge shuffling.
to construct the null distribution, we perform emdn on the randomized networks.
the empirical p-value of module is calculated as the probability of the random module having the observed score or greater by chance.
p-values are corrected for multiple testing using bh correction.
the cutoff of p-value is  <dig> , and it is statistically considered.
functional analysis of epigenetic modules
to assess the functional relevance of the epigenetic modules, we perform the gene ontology  enrichment analysis by using the hypergeometric test.
we obtain enrichment lists with bh corrected p-value .
features of svm on functional epigenetic modules
given a module c, we obtain the methylation level of each gene across all patient samples  <cit> , denoted by x
ij for the i-th gene and j-th patient.
for each sample j, the activity score of the k-th module is defined as the average gene methylation of all genes within the module, i.e.,
<dig> xc=∑i∈cxij/|c|.
where |c| is the number of genes in module c. a feature vector is constructed for each module.
simulated multiple networks
the simulated networks were generated following and the true classification of genes into clusters is known  <cit> .
we simulated three networks, each of which has  <dig> nodes.
in each network there are eight clusters of equal size.
the parameter controls the noise level of a network by controlling the ratio of intra-cluster edges to inter-cluster edges that are connected to a node.
the degree of each node is fixed to  <dig>  as mixed parameter increases from  <dig> to  <dig> , the detection of clusters in the networks becomes increasingly difficult.
the multiple networks contain two artificial networks, where the one network with noise level  <dig>  and the other with noise level from  <dig>  to  <dig> .
normalized mutual information
the normalized mutual information   <cit>  is based on the confusion matrix n whose rows correspond to the real modules in standard partition p
∗ and the columns correspond to the modules in obtained partition p. the element n
ij is the number of vertices overlapped by the i-th real and j-th predicted module.
the nmi is defined as
<dig> nmip,p∗=−2∑i=1|p|∑j=1|p∗|nijlognijnni.n.j∑i=1|p|ni.logni.n+∑i=1|p∗|n.jlogn.jn,
where |p| is the number of module in p and n
i. is the sum of the i-th row of matrix.
RESULTS
performance benchmarking of algorithms on artificial networks
for a comparative analysis, three algorithms are selected, including consensus clustering   <cit> , multiple-modularity method   <cit>  and spectral clustering   <cit> .
the spec method is not suitable for multiple networks.
therefore, the method is extended by using the consensus strategy  <cit> .
these algorithms are used in the artificial networks to test the performance.
figure 2
a is the heatmap of multiple networks where the common modules are surrounded by the dashed line.
we used the normalized mutual information  as a measure to quantify the performance.
before giving the performance of the compared algorithms, we first investigate how the parameter affects the performance of emdn.
the only parameter involved is the number of seeds.
the result is illustrated in fig.
2
b, where the number of seeds increased from  <dig> to 5%, and the accuracy increases dramatically.
the result indicates that 5% is the optimal value.
the accuracy of various algorithms on the artificial networks is shown in fig.
2
c, where emdn outperforms the others and the molti algorithm is better than the csc and spec methods.
fig.
<dig> performance of various algorithms on the artificial networks  heatmap of artificial multiple networks with noise level  <dig> ,  the number of seeds that affects the performance of emdn, and  performance of the compared algorithms on the artificial networks with various noise levels
fig.
<dig> the correlation between the methylation and expression of genes and modules is associated with breast cancer metastasis.
a heatmaps of genes and pathways between normal and tumor samples for dna methylation  and gene expression ,  indicator for correlation between the methylation and expression, where red denotes positively correlated, green for a negative correlation and white for noncorrelation, and  comparison between pairwise correlation coefficient among modules and random pathways based on methylation  and gene expression
negatively and positively correlated genes and pathways are common in breast cancer
prior to giving the performance of emdn algorithm, we first investigate the existence of positively correlated genes and pathways for breast cancer.
we select the acknowledged genes and pathways associated with breast cancer metastasis from ref.
<cit> , including  <dig> genes and  <dig> pathways as shown in fig.
3
a. the correlation between gene promoter methylation and gene expression is presented in fig.
3
b, where  <dig> genes are negatively correlated ,  <dig> positively correlated  and  <dig> uncorrelated .
the result is consistent with that presented in ref.
<cit> , where the positively correlated genes are also important in breast cancer metastasis.
moreover, the pathways can be categorized into two groups: pathways only with negatively correlated genes as well as pathways mixed with both positively and negatively correlated genes.
for example, the mek and growth factor receptor pathways are the first group, whereas the pi3k/akt/mtor and cell cycle are classified under the second group.
moreover, we find that the pi3k/akt/mtor and cell cycle are significantly positively correlated .
it is reasonable because majority of genes within pathway are positively correlated, for example, the coefficients are  <dig> , – <dig> , and  <dig>  for ccnd <dig>  cdkn2a, and cdk <dig>  respectively.
the result indicates the existence of positively correlated genes and pathways for breast cancer, which is also the main focus of this paper.
whether the positively and negatively correlated pathways are ubiquitous for breast cancer in various databases, including the kegg  <cit> , reactome  <cit> , biocart  <cit> .
we find that the pathways in all databases are both negatively and positively correlated although the majority of the pathways are negatively correlated .
specifically,  <dig> % of the pathways in biocart  are negatively correlated with coefficient ≤ – <dig>  , whereas  <dig> % of the pathways in biocart  are positively correlated with coefficient ≥  <dig> .
to check the significance of the percentage of correlated pathways, we compare the fold change of pathway percentage whose absolute value of correlation coefficient is equal or greater than  <dig>  with the size-matched random pathways .
the fold changes are  <dig>  ,  <dig>  , and  <dig>   for negatively correlated pathways, whereas the fold changes are  <dig>  ,  <dig>  ,  <dig>   for postively correlated pathways.
both positively and negatively correlated pathways are significantly higher than the random modules .
fig.
<dig> pathway expression and methylation correlation.
a the fold change of pathways in various databases in terms of correlation between gene expression and methylation, where ∗ denotes a significant difference when compared to the random modules , and  the hand <dig> modules, where node colours denote the differential dna methylation statistics as indicated.
border colors denote the differential expression statistics.
hexagon genes are recognized by emdn
the gene level correlation between methylation and expression profiles is also calculated.
as shown in additional file 2: figure s1a, the gene level methylation-expression correlations are both positive and negative.
for example,  <dig> genes with coefficients less than – <dig>  are observed, whereas only  <dig> genes with coefficients greater than  <dig>  are observed.
performance benchmarking of algorithms on hand <dig> module
to evaluate the performance of emdn algorithm, we compare it with several state-of-the-art, including consensus clustering   <cit> , multiple-modularity method   <cit> , fem  <cit> , epimod algorithm  <cit> .
because epimod is designed for single network, we applied it to differential comethylation and coexpression network respectively, denoted by epimod-meth and epimod-exp.
we employed the hand <dig> module as benchmark, which has  <dig> genes centered at hand <dig> gene .
the epimod-meth and epimod-exp algorithms cannot discover the module, while the rest algorithms discover the module.
the result demonstrates that the integrative analysis of methylation and expression data is promising for functional epigenetic modules.
the module recognized by emdn contains  <dig> genes, in which  <dig> genes are from hand <dig> module , whereas the accuracy of fem, csc and molti is  <dig> %,  <dig> % and  <dig> %.
the reason why the fem and emdn algorithms outperform csc and molti is that both css and molti are not designed for epigenetic modules.
the results indicates that multiple networks based strategy is more effective than the single network based approaches.
to fully characterize the performance, we compare fem and emdn on the simulated hand <dig> module  by using the accuracy , positive predictive value , false discovery rate  and false positive rate .
the comparison between fem and emdn is in table  <dig>  emdn is better than fem on acc, fpr, ppv as well as fdr.
two possible reasons are presented to explain why emdn is better than fem.
the first reason is that the multiple network model is a better way to characterize the functional epigenetic module than the single network based integration strategy because it avoids specifying the correlation between gene expression and methylation.
the second reason is that the module search strategy used in emdn is effective.
performance benchmarking of emdn on tcga breast cancer networks
next, we compare the algorithms using a compendium of  <dig> samples from tcga.
the differential comethylation  network has  <dig>  genes and  <dig>  edges.
we identify  <dig>   <dig>   <dig> and  <dig> modules using emdn, epimod-exp, epimod-meth and emdn, respectively.
no gold-standard exists for epigenetic modules in breast cancer.
thus, we compare the algorithms from three aspects.
first, the correlation coefficients between the methylation and expression of the modules are used to validate the performance because the ultimate goal of the algorithms is to discover methylated gene modules.
the higher the coefficient is, the more likely the methylated module is.
we compare the absolute value of correlation coefficients for epigenetic modules obtained each algorithm as shown in fig.
5
b. we conclude that the coefficients of the modules obtained by emdn are significantly higher than those obtained by other algorithms .
the result demonstrate that the proposed algorithm is more effective in discovering methylated gene modules.
fig.
<dig> performance benchmarking on tcga data.
a percentage of negatively and positively correlated modules, where green denotes a positive correlation and blue for negative correlation,  distribution of absolute values of module correlation between methylation and expression, where p-values are computed using the student t-test,  specificity of the methods.
gene modules found by each method are evaluated by a set of gold-standard pathway annotation.
specificity is defined as the fraction of gene modules significantly enriched for genes of some reference sets, and  sensitivity of the methods.
sensitivity is defined as the fraction of reference sets significantly enriched for genes of some modules.
enrichment p-values are computed using hypergeometric test and corrected by benjamin-hochberg, where * denotes p< <dig>
second, both positively and negatively correlated modules for breast cancer are shown.
moreover, whether all these algorithms can discover both negatively and positively correlated modules is questioned.
the percentages of positively and negatively correlated modules for each algorithm are shown in fig.
5
a. we conclude that only fem algorithm cannot recognize positively correlated modules.
the reason is that fem assumes the negative correlation between gene expression and methylation.
finally, we check how well the modules are enriched in the known pathways.
multiple reference pathway annotations are used, including go  <cit> , kegg  <cit> , reactome  <cit> , biocart  <cit> , and canonical pathways from msigdb  <cit> .
to evaluate the performance, we use specificity and sensitivity to quantify the accuracy, as shown in fig.
5
c and d. based on the results, we conclude that emdn achieves significantly higher specificity when evaluated using all reference sets while a comparable sensitivity is maintained .
epigenetic modules serve as biomarkers that improve the accuracy of breast cancer subtype prediction
the genes or pathways has been used to improve the prognosis of breast cancer  <cit> .
we hypothesize that the epigenetic modules may also likely be considered as biomarkers that predict luminal a/b, her <dig> and basal-like breast cancer subtypes.
we use the support vector machine  package for multiple classes classification  <cit> .
as a baseline comparison, we first compare the classification accuracy using the following feature sets: modules obtained by various algorithms and size-matched set of random modules.
svm is used for classification by using methylation profiles , where the area under curve  is employed to measure performance.
the results on tcga data using fivefold cross validation are shown in fig.
6
a, suggesting that the modules obtained by emdn are more discriminative than others on breast cancer subtype prediction.
the emdn had significantly higher auc .
fig.
<dig> epigenetic modules as biomarkers for tumor prediction.
a classification of tcga data.
roc curves for predicting breast cancer tumor using different feature sets, including random modules  validation of prediction on external data.
the classification model is trained on tcga data and tested on external data
to check the possibilities that the above result is biased because of parameter selection, we perform additional analysis by varying each of these parameters.
the result is consistent when we employ a 3/10-fold cross validation .
to check out the possibility that the confounding factors in tcga dataset contribute to the classification accuracy, we evaluate the performance of the svm classifiers  using an external data   <cit> .
the result indicates that the similar tendency is consistent with that in the tcga dataset as shown in fig.
6
b. the results indicate that the performance is not due to the hidden confounding factors in the tcga data .
finally, we investigate the possible reasons why emdn outperforms fem.
we quantify the importance of each epigenetic module, which is defined as the auc of svm based on the single module for subtype prediction.
we compare the feature importance between positively and negatively correlated modules.
there is no significant difference between the positively and negatively correlated modules .
both of positive and negative correlation epigenetic modules are discriminative for cancer tumor prediction.
epigenetic modules are associated with survival time of patients
genes and modules are associated with patient survival time in breast cancer  <cit> .
we hypothesize that the epigenetic modules are also associated with the clinical outcome using gene methylation.
the multivariable cox proportional hazard models are used to predict the survival time .
for each module, we predict the patient survival time by using the methylation  and the patients are segregated well into high- and low-risk groups according to the patient survival time.
we find that both the positively and negatively correlated modules divide the patients into two groups whose survival time is significant different .
for example, the module  <dig> and  <dig> are negatively correlated and the patients are divided into two groups with significantly different survival time .
the module  <dig> and  <dig> are positively correlated and the patients are divided into two groups with significantly different survival time .
in module  <dig>  the enriched go term is immune response  with genes fcer1a, cd1b, cd1a and ms4a <dig>  etc.
the results further demonstrated that both the positively and negatively correlated modules are critical for patient survival analysis.
fig.
<dig> kaplan-meier survival analysis for patients based on the epigenetic modules.
negatively correlated modules for survival time of patients:  module  <dig> with a p-value of  <dig>  obtained using the log-rank test,  module  <dig> with a p-value of  <dig> .
positively correlated modules for survival time of patients:  module  <dig> with a p-value of  <dig> ,  module  <dig> with a p-value of  <dig>
moreover, seven out of  <dig> modules obtained by emdn are significantly associated with the survival time of patients, whereas  <dig> out of  <dig> random modules with similar sizes are significantly associated with survival time of patients .
comparison between coexpression/comethylation networks and differential networks
it is necessary to validate the possibility to replace the differential networks with the coexpression and comethylation networks for emdn.
we identify  <dig> modules by using the emdn algorithm based on the coexpression and comethylation networks .
we first check the module methylation expression correlation coefficients between co-exp/meth modules and differential modules.
the differential modules have significantly higher coefficients than the co-exp/meth modules .
then, we check the sensitivity and specificity of the two groups of modules, showing that the differential modules are more enriched by the known pathways than the co-exp/meth modules .
finally, we compare the two groups of modules on breast cancer subtype prediction, and the result demonstrates that the differential modules are more discriminative than the co-exp/meth modules .
the result indicates that the differential network is preferred than coexpression  network.
CONCLUSIONS
exploring the functional epigenetic patterns is critical for understanding the mechanisms of biological processes.
recent technology has made it possible to perform simultaneously multi-platform genomic profiling of biological samples, including dna methylation and gene expression.
however, the systematic and integrative analysis of heterogeneous data for discovering biologically relevant patterns is currently scarce.
currently, all the available methods integrate the methylation and expression data into a scaffold network, such as protein interaction, which requires specifying the correlation between gene expression and methylation.
however, the specification of the correlation is unreasonable because there are both positive and negative correlation.
in this study, a novel strategy is presented to characterize functional epigenetic modules by using heterogeneous differential networks.
the functional epigenetic module discovery corresponds to find common modules in differential coexpression and comethylation networks.
overall, emdn possesses several unique advantages.
it provides a novel strategy for the integrative analysis of methylation and expression data.
it is more effective and accurate than modern methods because correlation specification is not required.
the method is easy to extend for other data, such as illumina 27k methylation, and chip-seq.
the basic concept of emdn should be modified in future work.
first, the algorithm is a generalized framework to any cohort of expression and methylation data, although this study uses breast cancer as a proof-of-principle.
second, data integration  might further enhance the identification of complicated molecular events involved in heterogeneous data.
additional files
additional file  <dig> supplementary figure s <dig>  distribution of pathway methylation and expression correlation coefficients for various databases:  negatively correlated pathways, and  positively correlated pathways.
additional file  <dig> supplementary figure s <dig>  positively and negatively correlated genes in breast cancer between methylation and expression:  distribution of positively and negatively correlated coefficients, and  pie chart for positively and negatively correlated genes.
additional file  <dig> supplementary figure s <dig>  breast cancer subtype prediction:  accuracy of compared algorithms using a 3-fold cross validation,  accuracy of compared algorithms using a 10-fold cross validation, and  feature importance of positively and negatively correlated modules in classifying the breast cancer subtypes.
additional file  <dig> supplementary figure s <dig>  comparison between coexpression/methylation and differential networks:  distributions of correlation coefficients of modules,  accuracy on classifying the breast cancer subtypes,  specificity of modules on the enrichment of known pathways, and  sensitivity of modules on the enrichment of known pathways.
abbreviations
hand2heart- and neural crest derivatives-expressed protein 2
keggkyoto encyclopedia of genes and genomes
svdsingular value decomposition
svmsupport vector machine.
tcga: the cancer genome atlas - cancer genome
tsstranscription start site

 dna methylation is a chemical modification of cytosine bases, which is critical for cellular differentiation, cell development and disease progression. for example, dna methylation directly inhibits the binding of transcription factors, and methylation aberrations either predispose to or result in disease progression. with biotechnological advancements, dna methylation is also considered a biomarker of epigenome analysis. multiple platforms, including reduced representation bisulfite sequencing, whole genome shotgun bisulfite sequencing, methylc-seq, and target capture bisulfite sequencing and dna methylation beadarray, have been developed to generate genome-wide dna methylation data. high-throughput technologies have generated large-scale genome-wide dna methylation profiles for various cancers and cell lines, providing great opportunities for revealing the epigenetic mechanisms. various approaches have been proposed on the basis of methylation profiles to extract dna methylation patterns. for example, fleischer et al. identified cpg probes associated with the survival time of breast cancer patients. hinoue et al. recognized four distinct subgroups in colorectal cancer by analyzing large-scale genome-wide dna methylation profiles. yang et al. discovered the comethylation modules across cell lines by using a weighted comethylation network. varley et al. identified dynamic dna methylation patterns across human cell lines. gevaert et al. designed methylmix to obtain pancancer clusters which reveal a novel similarity across various cancers. further details have also been provided. although dna methylation patterns have been extensively investigate, open questions have yet to be answered. for example, it is acknowledged that dna methylation aberrations cause diseases by mediating gene expression. however, how the dna methylation regulates gene expression remains unknown. epigenetic modules is critical for revealing the epigenetic mechanisms of cancers. thus, algorithms for functional epigenetic modules should be established by simultaneously analyzing methylation and gene expression data. however, the design of algorithms for the integrative analysis of dna methylation and gene expression is highly nontrivial because of two reasons. first, integrative analysis requires large-scale sample matched methylation and gene expression profiles. second, dna methylation patterns are difficult to determine because the relationship between methylation and expression is unknown. for instance, it is acknowledged that the correlation between promoter methylation and gene expression is negative. however, the recent evidence indicates that the positive correlation is also presented. regarding the first issue, the consortium, such as tcga, has generated sample matched dna methylation and gene expression data for various cancers. thus, functional epigenetic module can be identified. regarding the second issue, researchers developed many algorithms, such as the epimod algorithms, based on the gene comethylation network by using methylation data and searching modules in the network. however, these algorithms are limited because they are solely based on either methylation or gene expression. consequently, these algorithms fail to obtain epigenetic modules precisely. to address this problem, researchers developed the fem algorithm that integrates dna methylation and gene expression into a protein interaction network by assuming the anti-correlation between dna methylation and gene expression. the algorithm successfully identifies a novel epigenetically deregulated hotspot and methylated gene modules. this finding indicates the superiority of integrative analysis in terms of revealing methylated gene modules. current algorithms are limited by heterogeneous data integrated into a network, and thus require the correlation between dna methylation and gene expression. however, the correlation between heterogeneous data is unknown and difficult to determine. recently, multiple networks have been widely used to characterize the complex biological patterns. for example, ma et al. used multiple networks to characterize the dynamic modules. cantini et al. proposed a novel multi-network-based strategy to integrate different layers of genomic information and developed algorithm to identify cancer driver genes. didier et al. assessed aggregation, consensus and multiplex-modularity approaches to detect communities from multiple networks. these algorithms indicate that characterizing biological patterns in multiple networks is more accurate than those in a single network. therefore, we propose an alternative method that characterizes the epigenetic modules by using multiple networks. we first construct the differential comethylation and coexpression networks, and define the common modules within multiple networks as the epigenetic modules. the proposed strategy is advantageous because it avoids specifying the correlation between heterogeneous data and the accuracy in the identification of epigenetic modules is improved. we also develop the emdn algorithm by simultaneously analyzing multiple networks. experimental results demonstrate that the proposed algorithm is more accurate than other algorithms. we further demonstrate that epigenetic modules enriched by the known pathways, serve as biomarkers to predict breast cancer subtypes and survival time of patients. the proposed model and algorithm provide an effective way for the integrative analysis of dna methylation and gene expression. fig. flowchart of the proposed algorithm. a the common modules in both gene expression and methylation networks are defined as the functional epigenetic modules, emdn consists of two major components: differential network construction and functional epigenetic module discovery. the differential network construction consists of two steps: first, the coexpression network is constructed using the gene expression profiles, then the p-values of gene expression between tumor and non-tumor samples are incorporated into the coexpression network to develop the differential networks. finally, the algorithm discovers the common modules in multiple differential networks methods the emdn consists of two components: network construction and epigenetic module discovery. step 1: differential network construction the development of a gene differential coexpression network involves two steps: 1) the development of a binary coexpression network and 2) edge weight assignment based on differential gene expression between normal and tumor conditions. we assume that the gene expression data represents normalized estimates of gene expression intensity and are summarized at the gene level, which includes rna-seq or expression data generated by using illumina beadchip or affymatrix arrays. the binary gene coexpression network is developed where the edge weight on a pair of genes is the absolute pearson coefficient of the gene expression profiles. to remove the indirect correlation due to a third gene, w the first order partial pearson correlation coefficient is used. finally, edges whose weights are equal or greater than the pre-defined threshold δ are selected. using the limma package, the p-values of the gene expression difference between normal and tumor samples for each gene are obtained. the p-values are corrected by benjamini-hochberg. then, the gene differential coexpression network is developed by using the p-values. specifically, the weight w ij on edge in the differential network is defined as wij=logpi+logpj1/22maxllogpl1/2if|cor|≥δ,0if|cor|<δ, where p i is the p-value of the differential expression for gene i, and c o r is the value of pearson correlation between the i-th and j-th gene based on the expression profiles. the genes that are co-expressed and significantly differentially expressed are assumed to have assigned heavy weights. in this study, the δ is. this value leads to the maximal number of genes connected in all the networks. analogously, the differential comethylation network is developed using the gene methylation profiles. step 2: discovering the epigenetic modules to infer the epigenetic modules, emdn is based our recently developed m-module algorithm that is designed for identifying common modules in multiple molecular interaction networks. emdn consists of three steps: seed prioritization, module search by seed expansion and refinement of candidate modules. the seed prioritization ranks genes in multiple networks by using the topological feature of the gene in the networks. specifically, for each network g k with the weighted adjacency matrix w k=n×n, a function f k:v→r is developed, where f k denotes the importance of the i-th gene in the corresponding network. the function is defined as fk=αak′fk+y, where ak′fk denotes the topological importance of nodes and y is the vector for the prior information. the parameter α controls the relative weight of the topological importance and prior knowledge. ak′ is normalized adjacency matrix, i.e., ak′=dk−1/2wkdk−1/ with d k=d i a g. we use the following iteration-based algorithm to obtain f k: fk=αak′fk+y, where t denotes the iteration, and fk = no prior information is used. thus, we set y= the iteration is terminated if there is no change between f and f. usually, only iterations are observed. we calculate the z-score for f k. after ranking the genes in all networks, the ranks of genes are obtained in multiple networks, denoted by r=[f f …,f m]. the final gene ranking is obtained by using the z-score of genes in multiple networks, i.e., the row sum of r. starting with each seed, the module search step iteratively adds genes whose addition causes the maximum decrease in the graph entropy-based objective function until no decrease in the objective function is observed. given a module c, the entropy function for network g k is defined as hk=−pklogpk−log, where pk=∑i,j∈cwijk/∑i∈cwijk. the entropy function for c across all networks is defined as h=∑k=12hk/|c|. the module search procedure expands the module by using eq.. each time the gene whose addition leads to the maximum decrease of h is included in c until no gene can improve h. during the refinement step, candidate modules whose sizes are smaller than are removed. to merge overlapping function modules, two modules whose jaccard index is greater than are merged. data we download gene expression and dna methylation data for breast cancer from tcga database. specifically, samples with matched level- illumina 450k methylation data and hiseq rsem gene-normalized rna-seq data are obtained with tumor and normal samples. for methylation data, the β signal of the probe is used, which is calculated as the methylated signal divided by the sum of the methylated and unmethylated signal. for the rna-seq data, the reads per kilobase of exon model per million mapped reads is used. the clinical information is also obtained from tcga. in all these datasets, probes with more than 30% missing values are removed, and probes with less than 30% missing values are imputed using the r package pamr. gene methylation profiles to assign dna methylation to a given gene, we follow the strategy in the ref.. specifically, the average value of the probes mapping within bp of the transcription start site is assigned to the gene. if no probes mapped within bp of the tss, we use the average value of probes mapping to the 1st exon of the gene. if such probes are also not available, we use the average value of probes mapping within bp of the tss. eigengenes of modules the eigengene of a module is defined as the first principal component based on singular value decomposition. in details, the gene expression matrix of a given module is denoted by x= where the index i corresponds to the module genes and the index j corresponds to the samples. the singular value decomposition of x is denoted by x=udvt, where the columns of the matrices u and v are the left- and right-singular vectors, respectively. the first column of v is the module expression eigengene. similarly, we obtain the module methylation eigengene by the gene methylation profiles. correlation of module methylation and expression the correlation between module methylation and expression is defined as the pearson correlation coefficient between module methylation and methylation eigengenes. the correlation is significant when bh correction is used with the cutoff. if the sign of pearson coefficient is greater than a positive correlation is observed, and negative if otherwise. survival analysis we use the function coxph to implement the cox proportional hazard model to analyze the association of methylation profile of each epigenetic module with the patient survival. we use the prognostic index to generate high- and low-risk patient groups indexi=∑c=1kβcxci where k is the number of cancer-specific modules, β c is the regression coefficient of the cox proportional hazard model for the c-th module and x ci is the average methylation level of genes within the c-th module in the i-th patient. patients are grouped into high- and low-risk groups based on the median of prognostic index. the survival difference between these two groups of patients is obtained by using the kaplan-meier estimator and log-rank method. hand module and simulation as the true module, we select the hand module, since the biological and clinical significance of the driver gene, hand has been extensively validated. to fully assess the sensitivity and specificity of the proposed method, the simulation model of hand module is also adapted, in which it simulates statistics of differential methylation and differential expression on the protein interaction network. the model bootstrapped statistics for the member genes of this module to come from the top and lower 5% statistics quartiles, with the statistics of the rest of the network genes bootstrapped from the middle 90% portion. for each simulation run, the accuracy is recorded. statistical significance of modules the statistical significance of modules is based on the null score distribution of the random modules generated using randomized networks. each network is completely randomized times by degree-preserved edge shuffling. to construct the null distribution, we perform emdn on the randomized networks. the empirical p-value of module is calculated as the probability of the random module having the observed score or greater by chance. p-values are corrected for multiple testing using bh correction. the cutoff of p-value is, and it is statistically considered. functional analysis of epigenetic modules to assess the functional relevance of the epigenetic modules, we perform the gene ontology enrichment analysis by using the hypergeometric test. we obtain enrichment lists with bh corrected p-value. features of svm on functional epigenetic modules given a module c, we obtain the methylation level of each gene across all patient samples, denoted by x ij for the i-th gene and j-th patient. for each sample j, the activity score of the k-th module is defined as the average gene methylation of all genes within the module, i.e., xc=∑i∈cxij/|c|. where |c| is the number of genes in module c. a feature vector is constructed for each module. simulated multiple networks the simulated networks were generated following and the true classification of genes into clusters is known. we simulated three networks, each of which has nodes. in each network there are eight clusters of equal size. the parameter controls the noise level of a network by controlling the ratio of intra-cluster edges to inter-cluster edges that are connected to a node. the degree of each node is fixed to as mixed parameter increases from to, the detection of clusters in the networks becomes increasingly difficult. the multiple networks contain two artificial networks, where the one network with noise level and the other with noise level from to. normalized mutual information the normalized mutual information is based on the confusion matrix n whose rows correspond to the real modules in standard partition p ∗ and the columns correspond to the modules in obtained partition p. the element n ij is the number of vertices overlapped by the i-th real and j-th predicted module. the nmi is defined as nmip,p∗=−2∑i=1|p|∑j=1|p∗|nijlognijnni.n.j∑i=1|p|ni.logni.n+∑i=1|p∗|n.jlogn.jn, where |p| is the number of module in p and n i. is the sum of the i-th row of matrix.  performance benchmarking of algorithms on artificial networks for a comparative analysis, three algorithms are selected, including consensus clustering, multiple-modularity method and spectral clustering. the spec method is not suitable for multiple networks. therefore, the method is extended by using the consensus strategy. these algorithms are used in the artificial networks to test the performance. figure 2 a is the heatmap of multiple networks where the common modules are surrounded by the dashed line. we used the normalized mutual information as a measure to quantify the performance. before giving the performance of the compared algorithms, we first investigate how the parameter affects the performance of emdn. the only parameter involved is the number of seeds. the result is illustrated in fig. 2 b, where the number of seeds increased from to 5%, and the accuracy increases dramatically. the result indicates that 5% is the optimal value. the accuracy of various algorithms on the artificial networks is shown in fig. 2 c, where emdn outperforms the others and the molti algorithm is better than the csc and spec methods. fig. performance of various algorithms on the artificial networks heatmap of artificial multiple networks with noise level, the number of seeds that affects the performance of emdn, and performance of the compared algorithms on the artificial networks with various noise levels fig. the correlation between the methylation and expression of genes and modules is associated with breast cancer metastasis. a heatmaps of genes and pathways between normal and tumor samples for dna methylation and gene expression, indicator for correlation between the methylation and expression, where red denotes positively correlated, green for a negative correlation and white for noncorrelation, and comparison between pairwise correlation coefficient among modules and random pathways based on methylation and gene expression negatively and positively correlated genes and pathways are common in breast cancer prior to giving the performance of emdn algorithm, we first investigate the existence of positively correlated genes and pathways for breast cancer. we select the acknowledged genes and pathways associated with breast cancer metastasis from ref., including genes and pathways as shown in fig. 3 a. the correlation between gene promoter methylation and gene expression is presented in fig. 3 b, where genes are negatively correlated, positively correlated and uncorrelated. the result is consistent with that presented in ref., where the positively correlated genes are also important in breast cancer metastasis. moreover, the pathways can be categorized into two groups: pathways only with negatively correlated genes as well as pathways mixed with both positively and negatively correlated genes. for example, the mek and growth factor receptor pathways are the first group, whereas the pi3k/akt/mtor and cell cycle are classified under the second group. moreover, we find that the pi3k/akt/mtor and cell cycle are significantly positively correlated. it is reasonable because majority of genes within pathway are positively correlated, for example, the coefficients are, –, and for ccnd cdkn2a, and cdk respectively. the result indicates the existence of positively correlated genes and pathways for breast cancer, which is also the main focus of this paper. whether the positively and negatively correlated pathways are ubiquitous for breast cancer in various databases, including the kegg, reactome, biocart. we find that the pathways in all databases are both negatively and positively correlated although the majority of the pathways are negatively correlated. specifically, % of the pathways in biocart are negatively correlated with coefficient ≤ –, whereas % of the pathways in biocart are positively correlated with coefficient ≥. to check the significance of the percentage of correlated pathways, we compare the fold change of pathway percentage whose absolute value of correlation coefficient is equal or greater than with the size-matched random pathways. the fold changes are,, and for negatively correlated pathways, whereas the fold changes are,, for postively correlated pathways. both positively and negatively correlated pathways are significantly higher than the random modules. fig. pathway expression and methylation correlation. a the fold change of pathways in various databases in terms of correlation between gene expression and methylation, where ∗ denotes a significant difference when compared to the random modules, and the hand modules, where node colours denote the differential dna methylation statistics as indicated. border colors denote the differential expression statistics. hexagon genes are recognized by emdn the gene level correlation between methylation and expression profiles is also calculated. as shown in additional file 2: figure s1a, the gene level methylation-expression correlations are both positive and negative. for example, genes with coefficients less than – are observed, whereas only genes with coefficients greater than are observed. performance benchmarking of algorithms on hand module to evaluate the performance of emdn algorithm, we compare it with several state-of-the-art, including consensus clustering, multiple-modularity method, fem, epimod algorithm. because epimod is designed for single network, we applied it to differential comethylation and coexpression network respectively, denoted by epimod-meth and epimod-exp. we employed the hand module as benchmark, which has genes centered at hand gene. the epimod-meth and epimod-exp algorithms cannot discover the module, while the rest algorithms discover the module. the result demonstrates that the integrative analysis of methylation and expression data is promising for functional epigenetic modules. the module recognized by emdn contains genes, in which genes are from hand module, whereas the accuracy of fem, csc and molti is %, % and %. the reason why the fem and emdn algorithms outperform csc and molti is that both css and molti are not designed for epigenetic modules. the results indicates that multiple networks based strategy is more effective than the single network based approaches. to fully characterize the performance, we compare fem and emdn on the simulated hand module by using the accuracy, positive predictive value, false discovery rate and false positive rate. the comparison between fem and emdn is in table emdn is better than fem on acc, fpr, ppv as well as fdr. two possible reasons are presented to explain why emdn is better than fem. the first reason is that the multiple network model is a better way to characterize the functional epigenetic module than the single network based integration strategy because it avoids specifying the correlation between gene expression and methylation. the second reason is that the module search strategy used in emdn is effective. performance benchmarking of emdn on tcga breast cancer networks next, we compare the algorithms using a compendium of samples from tcga. the differential comethylation network has genes and edges. we identify and modules using emdn, epimod-exp, epimod-meth and emdn, respectively. no gold-standard exists for epigenetic modules in breast cancer. thus, we compare the algorithms from three aspects. first, the correlation coefficients between the methylation and expression of the modules are used to validate the performance because the ultimate goal of the algorithms is to discover methylated gene modules. the higher the coefficient is, the more likely the methylated module is. we compare the absolute value of correlation coefficients for epigenetic modules obtained each algorithm as shown in fig. 5 b. we conclude that the coefficients of the modules obtained by emdn are significantly higher than those obtained by other algorithms. the result demonstrate that the proposed algorithm is more effective in discovering methylated gene modules. fig. performance benchmarking on tcga data. a percentage of negatively and positively correlated modules, where green denotes a positive correlation and blue for negative correlation, distribution of absolute values of module correlation between methylation and expression, where p-values are computed using the student t-test, specificity of the methods. gene modules found by each method are evaluated by a set of gold-standard pathway annotation. specificity is defined as the fraction of gene modules significantly enriched for genes of some reference sets, and sensitivity of the methods. sensitivity is defined as the fraction of reference sets significantly enriched for genes of some modules. enrichment p-values are computed using hypergeometric test and corrected by benjamin-hochberg, where * denotes p< second, both positively and negatively correlated modules for breast cancer are shown. moreover, whether all these algorithms can discover both negatively and positively correlated modules is questioned. the percentages of positively and negatively correlated modules for each algorithm are shown in fig. 5 a. we conclude that only fem algorithm cannot recognize positively correlated modules. the reason is that fem assumes the negative correlation between gene expression and methylation. finally, we check how well the modules are enriched in the known pathways. multiple reference pathway annotations are used, including go, kegg, reactome, biocart, and canonical pathways from msigdb. to evaluate the performance, we use specificity and sensitivity to quantify the accuracy, as shown in fig. 5 c and d. based on the results, we conclude that emdn achieves significantly higher specificity when evaluated using all reference sets while a comparable sensitivity is maintained. epigenetic modules serve as biomarkers that improve the accuracy of breast cancer subtype prediction the genes or pathways has been used to improve the prognosis of breast cancer. we hypothesize that the epigenetic modules may also likely be considered as biomarkers that predict luminal a/b, her and basal-like breast cancer subtypes. we use the support vector machine package for multiple classes classification. as a baseline comparison, we first compare the classification accuracy using the following feature sets: modules obtained by various algorithms and size-matched set of random modules. svm is used for classification by using methylation profiles, where the area under curve is employed to measure performance. the results on tcga data using fivefold cross validation are shown in fig. 6 a, suggesting that the modules obtained by emdn are more discriminative than others on breast cancer subtype prediction. the emdn had significantly higher auc. fig. epigenetic modules as biomarkers for tumor prediction. a classification of tcga data. roc curves for predicting breast cancer tumor using different feature sets, including random modules validation of prediction on external data. the classification model is trained on tcga data and tested on external data to check the possibilities that the above result is biased because of parameter selection, we perform additional analysis by varying each of these parameters. the result is consistent when we employ a 3/10-fold cross validation. to check out the possibility that the confounding factors in tcga dataset contribute to the classification accuracy, we evaluate the performance of the svm classifiers using an external data. the result indicates that the similar tendency is consistent with that in the tcga dataset as shown in fig. 6 b. the results indicate that the performance is not due to the hidden confounding factors in the tcga data. finally, we investigate the possible reasons why emdn outperforms fem. we quantify the importance of each epigenetic module, which is defined as the auc of svm based on the single module for subtype prediction. we compare the feature importance between positively and negatively correlated modules. there is no significant difference between the positively and negatively correlated modules. both of positive and negative correlation epigenetic modules are discriminative for cancer tumor prediction. epigenetic modules are associated with survival time of patients genes and modules are associated with patient survival time in breast cancer. we hypothesize that the epigenetic modules are also associated with the clinical outcome using gene methylation. the multivariable cox proportional hazard models are used to predict the survival time. for each module, we predict the patient survival time by using the methylation and the patients are segregated well into high- and low-risk groups according to the patient survival time. we find that both the positively and negatively correlated modules divide the patients into two groups whose survival time is significant different. for example, the module and are negatively correlated and the patients are divided into two groups with significantly different survival time. the module and are positively correlated and the patients are divided into two groups with significantly different survival time. in module the enriched go term is immune response with genes fcer1a, cd1b, cd1a and ms4a etc. the results further demonstrated that both the positively and negatively correlated modules are critical for patient survival analysis. fig. kaplan-meier survival analysis for patients based on the epigenetic modules. negatively correlated modules for survival time of patients: module with a p-value of obtained using the log-rank test, module with a p-value of. positively correlated modules for survival time of patients: module with a p-value of, module with a p-value of moreover, seven out of modules obtained by emdn are significantly associated with the survival time of patients, whereas out of random modules with similar sizes are significantly associated with survival time of patients. comparison between coexpression/comethylation networks and differential networks it is necessary to validate the possibility to replace the differential networks with the coexpression and comethylation networks for emdn. we identify modules by using the emdn algorithm based on the coexpression and comethylation networks. we first check the module methylation expression correlation coefficients between co-exp/meth modules and differential modules. the differential modules have significantly higher coefficients than the co-exp/meth modules. then, we check the sensitivity and specificity of the two groups of modules, showing that the differential modules are more enriched by the known pathways than the co-exp/meth modules. finally, we compare the two groups of modules on breast cancer subtype prediction, and the result demonstrates that the differential modules are more discriminative than the co-exp/meth modules. the result indicates that the differential network is preferred than coexpression network.  exploring the functional epigenetic patterns is critical for understanding the mechanisms of biological processes. recent technology has made it possible to perform simultaneously multi-platform genomic profiling of biological samples, including dna methylation and gene expression. however, the systematic and integrative analysis of heterogeneous data for discovering biologically relevant patterns is currently scarce. currently, all the available methods integrate the methylation and expression data into a scaffold network, such as protein interaction, which requires specifying the correlation between gene expression and methylation. however, the specification of the correlation is unreasonable because there are both positive and negative correlation. in this study, a novel strategy is presented to characterize functional epigenetic modules by using heterogeneous differential networks. the functional epigenetic module discovery corresponds to find common modules in differential coexpression and comethylation networks. overall, emdn possesses several unique advantages. it provides a novel strategy for the integrative analysis of methylation and expression data. it is more effective and accurate than modern methods because correlation specification is not required. the method is easy to extend for other data, such as illumina 27k methylation, and chip-seq. the basic concept of emdn should be modified in future work. first, the algorithm is a generalized framework to any cohort of expression and methylation data, although this study uses breast cancer as a proof-of-principle. second, data integration might further enhance the identification of complicated molecular events involved in heterogeneous data. additional files additional file supplementary figure s distribution of pathway methylation and expression correlation coefficients for various databases: negatively correlated pathways, and positively correlated pathways. additional file supplementary figure s positively and negatively correlated genes in breast cancer between methylation and expression: distribution of positively and negatively correlated coefficients, and pie chart for positively and negatively correlated genes. additional file supplementary figure s breast cancer subtype prediction: accuracy of compared algorithms using a 3-fold cross validation, accuracy of compared algorithms using a 10-fold cross validation, and feature importance of positively and negatively correlated modules in classifying the breast cancer subtypes. additional file supplementary figure s comparison between coexpression/methylation and differential networks: distributions of correlation coefficients of modules, accuracy on classifying the breast cancer subtypes, specificity of modules on the enrichment of known pathways, and sensitivity of modules on the enrichment of known pathways. abbreviations hand2heart- and neural crest derivatives-expressed protein 2 keggkyoto encyclopedia of genes and genomes svdsingular value decomposition svmsupport vector machine. tcga: the cancer genome atlas - cancer genome tsstranscription start site 
to address this issue, we present a novel multiple network framework for epigenetic modules, namely, epigenetic module based on differential networks algorithm, by simultaneously analyzing dna methylation and gene expression data. epigenetic modules also estimate the survival time of patients, and this factor is critical for cancer therapy. methylationnetwork biologymultiple networksepigenetic moduleissue-copyright-statement© the author 2017 with the increase in the amount of dna methylation and gene expression data, the epigenetic mechanisms of cancers can be extensively investigate. available methods integrate the dna methylation and gene expression data into a network by specifying the anti-correlation between them. the proposed model and algorithm provide an effective method for the integrative analysis of dna methylation and gene expression. these modules can serve as bio-markers to predict breast cancer subtypes by using methylation profiles, where positively and negatively correlated modules are of equal importance in the classification of cancer subtypes. the emdn algorithm prevents the specification of the correlation between methylation and expression. the accuracy of emdn algorithm is more efficient than that of modern approaches. on the basis of the cancer genome atlas breast cancer data, we observe that the emdn algorithm can recognize positively and negatively correlated modules and these modules are significantly more enriched in the known pathways than those obtained by other algorithms.
 with the increase in the amount of dna methylation and gene expression data, the epigenetic mechanisms of cancers can be extensively investigate. available methods integrate the dna methylation and gene expression data into a network by specifying the anti-correlation between them. however, the correlation between methylation and expression is usually unknown and difficult to determine.  to address this issue, we present a novel multiple network framework for epigenetic modules, namely, epigenetic module based on differential networks algorithm, by simultaneously analyzing dna methylation and gene expression data. the emdn algorithm prevents the specification of the correlation between methylation and expression. the accuracy of emdn algorithm is more efficient than that of modern approaches. on the basis of the cancer genome atlas breast cancer data, we observe that the emdn algorithm can recognize positively and negatively correlated modules and these modules are significantly more enriched in the known pathways than those obtained by other algorithms. these modules can serve as bio-markers to predict breast cancer subtypes by using methylation profiles, where positively and negatively correlated modules are of equal importance in the classification of cancer subtypes. epigenetic modules also estimate the survival time of patients, and this factor is critical for cancer therapy.  the proposed model and algorithm provide an effective method for the integrative analysis of dna methylation and gene expression. the algorithm is freely available as an r-package at https://github.com/william0701/emdn. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords methylationnetwork biologymultiple networksepigenetic moduleissue-copyright-statement© the author 2017 