NEpiC: a network-assisted algorithm for epigenetic studies using mean and variance combined signals

Abstract

DNA methylation plays an important role in many biological processes. Existing epigenome-wide association studies (EWAS) have successfully identified aberrantly methylated genes in many diseases and disorders with most studies focusing on analysing methylation sites one at a time. Incorporating prior biological information such as biological networks has been proven to be powerful in identifying disease-associated genes in both gene expression studies and genome-wide association studies (GWAS) but has been under studied in EWAS. Although recent studies have noticed that there are differences in methylation variation in different groups, only a few existing methods consider variance signals in DNA methylation studies. Here, we present a network-assisted algorithm, NEpiC, that combines both mean and variance signals in searching for differentially methylated sub-networks using the protein–protein interaction (PPI) network. In simulation studies, we demonstrate the power gain from using both the prior biological information and variance signals compared to using either of the two or neither information. Applications to several DNA methylation datasets from the Cancer Genome Atlas (TCGA) project and DNA methylation data on hepatocellular carcinoma (HCC) from the Columbia University Medical Center (CUMC) suggest that the proposed NEpiC algorithm identifies more cancer-related genes and generates better replication results.

INTRODUCTION

DNA methylation plays critical roles in many biological activities, especially in the carcinogenesis process. Two common kinds of aberrant methylation in cancers are regional hypermethylation and global hypomethylation. Hypermethylation within promoter regions of genes, which may lead to the silence of associated genes, is known to be associated with various kinds of cancers, such as liver (1), renal (2), colorectal (3) and endometrial (4) cancers. Global hypomethylation mainly affects intergenic regions of the genome and may increase chromosomal instability (5).

Many epigenome-wide association studies (EWAS) have successfully identified aberrantly methylated genes in cancers (6–8) with most studies focusing on analyzing DNA methylation sites one at a time. Several methods have being developed that consider correlations among sites on a gene or correlations among genes in a pathway (9,10) using penalized regression models. In genome-wide association studies (GWAS) or gene expression studies, incorporating prior biological information has been proven to be a more effective way to identify disease-associated single nucleotide polymorphisms (SNPs) or genes that are enriched with stronger association signals and higher biological relevance (11–15). Those methods prioritize candidate SNPs or genes by incorporating prior biological information such as gene annotations, biological pathways or protein–protein interaction (PPI) networks. More specifically, network-assisted methods overlay genetic or gene expression signals onto a biological network and search for sub-networks (modules) with GWAS data or gene expression data. Jia et al. (13) proposed a dense module searching method for GWAS (dmGWAS) that searches for modules that are enriched with genes of higher significances (low P-values) within a PPI network and showed that dmGWAS is more powerful in identifying disease related genes than other methods that do not incorporate network information. There are also some EWAS studies that incorporate biological network information (16,17). Another feature of DNA methylation measures that was recently observed is the higher variation in cancer tissues than in normal tissues across human cancer types (18). A few methods for DNA methylation data that consider differences in variances between two experimental conditions have already been developed (19–21). However, there is no method that incorporates both the prior biological information such as the network information and variance signals in DNA methylation studies. In this paper, we propose the NEpiC algorithm, a Network-assisted algorithm for Epigenetic studies that uses mean and variance Combined signals in searching for differentially methylated sub-networks in a PPI network.

In the proposed NEpiC algorithm, we first compute mean and variance signal scores for a CpG site and then summarize the two scores with weights to create a combined score for the CpG site. We then extract the gene-level score out of all CpG sites on a gene. Finally, using a PPI network, we search for dense modules that are enriched for genes with large gene-level scores with a greedy search algorithm. We conducted simulation studies to show the performance of the proposed NEpiC algorithm compared to methods that either do not use the biological network information or do not use variance signals in searching for differentially methylated genes. We applied NEpiC to the 450K DNA methylation datasets of tumor and adjacent normal tissues of breast invasive carcinoma (BRCA) and liver hepatocellular carcinoma (LIHC) from the Cancer Genome Atlas (TCGA) project as well as an independent 450K DNA methylation data of tumor and adjacent normal tissues of hepatocellular carcinoma (HCC) from the Columbia University Medical Center (CUMC) (8). The results show that the proposed NEpiC algorithm which uses the biological network information among genes and both mean and variance signals at each CpG site identifies more cancer-related genes and generates better replication results than methods that do not consider both pieces of information.

MATERIALS AND METHODS

Since matched case–control designs with tumor and adjacent normal tissues are frequently used in DNA methylation studies of cancer, we focused on studies with a matched case–control design here. However, the proposed NEpiC algorithm is readily modified and applied to other types of designs. There are three steps in the proposed NEpiC algorithm: (i) constructing site-level and gene-level signal scores using DNA methylation data; (ii) searching modules on the PPI network with the guide of gene-level scores; (iii) prioritizing modules and candidate genes; (iv) validating identified modules using permutations. Figure 1 displays the pipeline of the proposed NEpiC algorithm.

Figure 1.

NEpiC algorithm pipeline.

Step 1: Constructing scores

Site-level mean and variance signal scores

We use the two-sided paired t-test and the one-sided Morgan–Pitman Test (22,23) to calculate P-values to test if the mean methylation measures are the same between tumor and adjacent normal groups and if the variance of the methylation measures in the tumor groups is greater than that in the adjacent normal group at CpG site i, which are denoted as |${p_{mi}}$| and|$\ {p_{vi}}$|⁠, respectively. The mean and variance signal scores at CpG site i are then defined as |${m_i} = {\Phi ^{ - 1}}( {1 - {p_{mi}}} )$|⁠, and |${v_i} = {\Phi ^{ - 1}}( {1 - {p_{vi}}} )$|⁠, where |$\Phi$| is the standard normal distribution function. We set those mean and variance scores that are smaller than zero (i.e., those sites with mean and variance P-values > 0.5) to be zero and delete the CpG sites whose mean and variance scores are both zero.

Site-level combined signal scores

Due to potentially different scales for site-level mean and variance signal scores, we define a combined signal score |${c_i}$|at CpG site i weighted by λ: |${c_i} = \lambda {m_i} + ( {1 - \lambda } ){v_i}$| to balance the mean and variance scores. We calculate λ as follows: we first define a ratio |${r_i} = \frac{{{v_i}}}{{{m_i} + {v_i}}}$|at each CpG site i; we then average across all sites on a gene to obtain the gene-level ratio; lastly, we average gene-level ratios across all genes genome-wide to obtain λ.

Gene-level combined signal scores

We choose the CpG site with the largest combined signal score |${c_i}$| to represent gene j, and denote this gene-level score with|${\rm{\ }}{g_j}$|⁠,|${\rm{\ }}j = 1, \ldots ,J$|⁠, where |$J$| is the total number of genes.

Step 2: Searching for modules

To search for modules enriched with genes of high gene-level signal scores, we define module scores S as follows:

\begin{equation*}S = \frac{{\mathop \sum \nolimits_{j = 1}^m {g_j}}}{{\sqrt m }},\end{equation*}

where m is the number of genes in a module. We use a greedy search algorithm with the following steps to search for modules. (i) Set a gene on the PPI network as the seed gene, which is considered as the starting module, and calculate the module score S. (ii) Identify the gene with the largest gene-level score |${g_j}$|from all genes that are the first order neighbors of the seed gene on the PPI network and add this gene to the starting module only if the module score increases by a predefined cutoff, e.g., >10% (13). (iii) Continue with the first-order genes of the genes in the current module and keep adding genes if the module score increases by >10%, otherwise stop the algorithm and save the current module as the module corresponding to the seed gene. (iv) Choose another gene on the PPI network as the seed gene and repeat the above steps until all genes on the PPI network have been considered as a seed gene.

Step 3: Prioritizing modules and candidate genes

We now have the same number of modules as the number of genes on the PPI network. We exclude the small modules of size smaller than five genes (13). Since modules with more genes have larger module scores S in general, to make modules of different sizes comparable, we normalize module scores according to their sizes through a resampling method. More specifically, for each module obtained, we randomly generate 100 000 modules of the same size from the PPI network and calculate their module scores. We then normalize the module score S and define a normalized module score |${S_N}$| as |${S_N} = \frac{{S - \mu }}{\sigma }$|⁠, where |$\mu$| and |$\sigma$| are the mean and standard deviation of those 100 000 module scores. We then rank the modules by their normalized module scores and select the top ranked modules. Genes that appear in the selected top-ranked modules are candidate genes that are potentially associated with the tumor status. We also rank the candidate genes using the frequencies that each candidate gene is selected by the selected top-ranked modules. We define those candidate genes that are selected by more than one module as the prioritized candidate genes.

Step 4: Validating identified modules using permutations

To control for the potential bias brought by different CpG site density or different gene sizes, we perform a permutation test on identified modules where we shuffle the tumor/adjacent normal labels within tumor/adjacent normal pairs and repeat the permutation procedures 100 times. We then calculate the permutated module scores |${S_{perm}}$| for the identified top-ranked modules. The permutation P-values of the identified top-ranked modules are calculated as follows:

\begin{equation*}{P_{perm}} = \frac{{\# \{ {S_{perm}} >S\} }}{{\# \left\{ {total\ permutation} \right\}}},\end{equation*}

where |$S$| is the observed module score.

Simulation studies

We conducted simulation studies to investigate the performance of the proposed NEpiC algorithm that uses the biological network information and both mean and variance signals in DNA methylation data, and compared the performance of NEpiC with that of the following three methods: (i) NEpi method that uses only mean signals of DNA methylation data and also incorporates the biological network information in searching for dense modules, which is an extension of dmGWAS (13); (ii) EpiC method that uses combined signals in both mean and variance of DNA methylation data but not the biological network information; and (iii) Epi method that uses only mean signals of DNA methylation data and not the biological network information.

Simulation settings

PPI network

We used the PPI network from the Protein Interaction Network Analysis platform (PINA) (24) and only maintained edges established with experimental evidence. We also excluded UBC gene from the PPI network which has a degree of connection of 9112, far greater than the rest of the genes on the network. The final PINA PPI network has 13 932 genes and 131 778 edges. We further trimmed the PPI network to keep genes that are also in the CUMC HCC methylation data after quality control steps (see Real Data Applications). We ended up with 12 630 genes and 116 772 edges.

Simulating site-level signal scores

We simulated site-level mean and variance signal scores based on real data from the CUMC HCC study where there are 195 259 CpG sites from 12 630 genes for 66 matched tumor and adjacent normal pairs after quality control steps which overlap with the genes in PINA PPI network.

To assign genes out of the total 12 630 genes on the PPI network as outcome-associated genes, we first chose the 25 genes that were considered as driver genes in a recent review paper on HCC (25) as the outcome-associated genes. We then randomly selected genes that are the first-order neighbors of the 25 driver genes as the outcome-associated genes with a probability 0.02. We considered the 25 driver genes plus the selected first-order outcome-associated genes as the ‘seed associated genes’, and randomly selected genes that are the first-order neighbors of those ‘seed associated genes’ as outcome-associated genes with a probability 0.02. We chose a probability of 0.02 such that after repeating this procedure 5–6 times we will have around 500 outcome-associated genes. This is because there are currently 572 mutated genes that have been causally implicated with cancers according to the Cancer Gene Census category (CGC, as of December 2015) (26). The rest around 12 130 genes are then set as not associated with the outcome.

To assign effect sizes to the outcome-associated genes simulated above, we first defined three groups of genes with different effect sizes, genes with strong, median, and weak effects, where the effect sizes are determined based on the results from the CUMC HCC data. Within each gene category, the mean and variance signal scores were then simulated with a bivariate normal distribution with means, variances and correlations mimicking results from the CUMC HCC data. We set the 25 driver genes to have strong effects. We randomly set half of the selected outcome-associated genes excluding the 25 driver genes to have median effects and the other half to have weak effects. Within each effect size group, we divided them into three subgroups. Using the gene group with a strong effect as the example, the three subgroups are: genes whose mean and variance signal scores are correlated; genes whose mean and variance signal scores are independent when mean signal scores have strong effects and variance signal scores have no effects; and genes whose mean and variance signal scores are independent when mean signal scores have no effects and variance signal scores have strong effects. The ratio of number of genes in these three subgroups is 2:4:1. See supplementary materials for more details of the simulation settings.

RESULTS

Simulation results

We simulated 100 datasets and applied the proposed NEpiC algorithm and the three compared methods, NEpi, EpiC, and Epi. For NEpiC and NEpi, we chose all genes in the top 1% of modules as candidates. For EpiC and Epi methods that do not use network information, we chose the same number of top-ranked genes as that in the top 1% of modules selected by NEpiC based on gene-level combined signal scores (EpiC method) or gene-level mean signal scores (Epi method). Part I in Table 1 shows the average number of candidate genes identified and the average number (percentage) of truly associated genes out of the candidate genes identified using the four methods and the enrichment P-value of the truly associated genes among candidate genes identified. The enrichment P-value was calculated using a hypergeometric distribution hyper(q, M, N, k), where q is the number of truly associated genes among identified candidate genes, M is the number of truly associated genes, N is the number of true null genes and k is the number of identified candidate genes. It shows that the proposed NEpiC algorithm identifies more truly associated genes with a higher percentage of truly associated genes out of the candidate genes identified and achieves a more significant enrichment P-value than the other three methods. The fact that the NEpiC algorithm outperforms the NEpi algorithm that uses the network information but not the variance signal in DNA methylation data and the EpiC algorithm that uses the variance signal in DNA methylation data but not the network information indicates the benefit of incorporating variance signals and the benefit of incorporating the network information, respectively. Part I of Table 1 also displays the average numbers of truly associated genes identified broken down by effect size categories. It is clear that the proposed NEpiC algorithm is more powerful in identifying genes with median or weak effects as expected.

Average numbers of candidate genes (Part I) and prioritized candidate genes (Part II) identified and truly associated genes among the candidate genes (Part I) and prioritized candidate genes (Part II) identified by NEpiC, NEpi, EpiC and Epi algorithms and the enrichment P-values based on 100 simulated datasets

Table 1.

Average numbers of candidate genes (Part I) and prioritized candidate genes (Part II) identified and truly associated genes among the candidate genes (Part I) and prioritized candidate genes (Part II) identified by NEpiC, NEpi, EpiC and Epi algorithms and the enrichment P-values based on 100 simulated datasets

	NEpiC	NEpi	EpiC	Epi
Part I
k = Number of candidate genes identified	151	149	151	151
q = Number of candidate genes that are truly associated (combining strong, median or weak effects) (%1)	21 (13.9%)	15 (10.1%)	15 (9.9%)	11 (7.3%)
Number of candidate genes that are truly associated with strong effect	11	9	14	11
Number of candidate genes that are truly associated with median effect	7	4	1	0
Number of candidate genes that are truly associated with weak effect	3	2	0	0
Enrichment P-value of truly associated genes among candidate genes2	|$1.37 \times {10^{ - 8}}$|	|$6.76 \times {10^{ - 5}}$|	|$7.94 \times {10^{ - 5}}$|	|$6.89 \times {10^{ - 3}}$|
Part II
k = Number of prioritized candidate genes identified	31	30	31	31
q = Number of prioritized candidate genes that are truly associated (combining strong, median or weak effects) (%3)	10 (32.3%)	8 (26.7%)	5 (16.1%)	4 (12.9)
Number of prioritized candidate genes that are truly associated with strong effect	9	7	5	4
Number of prioritized candidate genes that are truly associated with median effect	1	1	0	0
Number of prioritized candidate genes that are truly associated with weak effect	0	0	0	0
Enrichment P-value of truly associated genes among prioritized candidate genes4	|$3.92 \times {10^{ - 9}}$|	|$4.87 \times {10^{ - 6}}$|	|$1.22 \times {10^{ - 3}}$|	|$4.15 \times {10^{ - 3}}$|

¹Percent truly associated genes out of the candidate genes identified.

²Enrichment P-value was calculated using the hypergeometric distribution: hyper(q, M, N, k), where q is the average number of truly associated genes among candidate genes identified, M is the average number of simulated truly associated genes (M = 443), N is the average number of simulated truly non-associated genes (N = 12 187), and k is the average number of candidate genes in the top 1% of modules.

³Percent truly associated genes out of the prioritized candidate genes identified.

⁴Enrichment P-value was calculated using the hypergeometric distribution: hyper(q, M, N, k), where q is the number of truly associated genes among prioritized candidate genes identified, M is the average number of simulated truly associated genes (M = 456), N is the average number of simulated truly non-associated genes (N = 12 174), and k is the average number of prioritized candidate genes in the top 1% of modules.

We further chose the candidate genes that were selected by more than one module among the top 1% of modules as the prioritized candidate genes using the NEpiC or NEpi algorithm. For the EpiC and Epi algorithms, we then chose the number of top ranked candidate genes as the same number of prioritized candidate genes identified using the proposed NEpiC algorithm (Table 1 Part II). Part II of Table 1 also displays the average number (percentage) of truly associated genes out of the prioritized candidate genes identified. There are 31 prioritized candidate genes on average using the proposed NEpiC algorithm (there are 30 prioritized candidate genes on average using the NEpi algorithm). The enrichment P-values of the truly associated genes among the prioritized candidate genes are more significant than that among the candidate genes for both the NEpiC and NEpi algorithms (Table 1). Moreover, the percentage of truly associated genes out of the prioritized candidate genes is higher than that out of the candidate genes for both the NEpiC and NEpi algorithms (Table 1). These suggest that the prioritization procedure using the selection frequencies might be a useful step to further improve the performance of the proposed NEpiC algorithm to identify truly outcome-associated genes.

Real data applications

We applied the proposed NEpiC algorithm to the TCGA BRCA, TCGA LIHC and CUMC HCC DNA methylation datasets. We removed probes with missing data in more than 70% of the samples and probes on the sex chromosomes and probes that are known SNP sites. We also removed probes with no gene annotations and required CpG coverage to be at least 95% per sample. Finally, we used Bioconductor package wateRmelon to correct for the type II probe bias (27).

TCGA BRCA data

After quality control steps, there are 229 655 CpG sites from 19 270 genes for 90 matched tumor and adjacent normal pairs in the TCGA BRCA data. Of those, 12 561 genes are also in the PPI network, which contains 115 928 edges. We then applied the proposed NEpiC algorithm using the PINA PPI network to the TCGA BRCA DNA methylation data.

Figure 2 shows the module scores of modules of different sizes before and after normalization using the proposed NEpiC algorithm. The module scores before normalization increase with module sizes as expected while the module scores after normalization are comparable. With the normalized module scores, we chose genes in the top 1% of modules as the candidate genes using NEpiC and NEpi, where there are 227 and 161 genes, respectively (Table 2 Part I). All the top 1% of modules identified by NEpiC and NEpi with the TCGA BRCA data have permutation P-values smaller than 0.0005. We then chose the top 227 genes with the strongest combined or mean signals using EpiC and Epi. Among the candidate genes identified by NEpiC, NEpi, EpiC and Epi, there are 16, 2, 11 and 1 genes that have been reported to be differentially methylated in cancers according to a cancer methylation database developed by combining text-mining and expert annotation (Pubmeth) (28) (Table 2 Part I). According to CGC (26), there are 26, 12, 10 and 11 genes that have been causally implicated with cancers, respectively (Table 2 Part I). NEpiC has the highest percentages of reported differentially methylated genes in cancers based on Pubmeth and causally implicated cancer genes based on CGC out of the candidate genes identified. This suggests that the proposed NEpiC algorithm that uses both the biological network information and the mean and variance signals is a more powerful method in identifying potentially cancer-related genes than the methods that ignore either the biological network information or variance signals in DNA methylation data. Between the list of 16 genes reported to be differentially methylated in cancers based on Pubmeth and the list of 26 genes reported to be causally implicated in cancers based on CGC, there are four genes in common. Although the number of genes in common is not large, we found that among the 22 causally implicated cancer genes identified by the CGC database only, nine genes were also reported to be aberrantly methylated in different cancers (5 in breast cancer and four in other cancers, all were published after the Pubmeth database was generated), seven genes were reported to be aberrantly expressed or mutated in breast cancer, and the remaining 6 genes were also reported to be aberrantly expressed or mutated in other cancers (Additional Table 1 for BRCA, Part I). We then conducted pathway enrichment analyses using WebGestalt (29), where we used the enrichment P-value of the ‘pathway in cancer’ from KEGG as the benchmark (Table 2 Part I). It shows that the proposed NEpiC algorithm achieves the most significant enrichment P-value while NEpi achieves the second most significant enrichment P-value.

Figure 2.

Module scores before and after normalization using the proposed NEpiC algorithm with the TCGA BRCA data.

Number of candidate genes (Part I) and prioritized candidate genes (Part II) identified in the TCGA BRCA data and number of reported differentially methylated genes in cancers and causally implicated cancer genes out of the candidate genes (Part I) and out of the prioritized candidate genes (Part II) identified according to PubMeth¹ and CGC² and enrichment P-value of ‘pathways in cancer’ from KEGG among the candidate genes (Part I) and the prioritized candidate genes (Part II) identified

Table 2.

Number of candidate genes (Part I) and prioritized candidate genes (Part II) identified in the TCGA BRCA data and number of reported differentially methylated genes in cancers and causally implicated cancer genes out of the candidate genes (Part I) and out of the prioritized candidate genes (Part II) identified according to PubMeth¹ and CGC² and enrichment P-value of ‘pathways in cancer’ from KEGG among the candidate genes (Part I) and the prioritized candidate genes (Part II) identified

	NEpiC	NEpi	EpiC	Epi
Part I
Number of candidate genes identified	227	161	227	227
Number (percentage) of reported differentially methylated genes in cancers out of the candidate genes according to Pubmeth (%3)	16 (7.0%)	2 (1.2%)	11 (4.8%)	1 (0.4%)
Number (percentage) of causal implicated cancer out of the candidate genes according to CGC (%4)	26 (11.5%)	12 (7.5%)	10 (4.4%)	11 (4.8%)
Enrichment P-value of ‘pathway in cancer’ among the candidate genes identified	|$1.96 \times {10^{ - 13}}$|	|$2.00 \times {10^{ - 4}}$|	NS5	|$5.00 \times {10^{ - 3}}$|
Part II
Number of prioritized candidate genes identified	68	27	68	68
Number (percentage) of reported differentially methylated genes in cancers out of the prioritized candidate genes according to Pubmeth (%6)	8 (11.8%)	1 (3.7%)	2 (2.9%)	1 (1.5%)
Number (percentage) of causally implicated cancer genes out of the prioritized candidate genes according to CGC (%7)	14 (20.6%)	3 (11.1%)	2 (2.9%)	3 (4.4%)
Enrichment P-value of ‘pathway in cancer’ among the prioritized candidate genes identified	|$1.42 \times {10^{ - 13}}$|	|$4.90 \times {10^{ - 5}}$|	NS	NS

¹There are 292 reported differentially methylated genes in cancers according to Pubmeth (28).

²There are 572 mutated genes that have been causally implicated with cancers according to the Cancer Gene Census category (CGC, as of December 2015) (26).

³Percent Pubmeth genes out of the candidate genes identified.

⁴Percent CGC genes out of the candidate genes identified.

⁵NS stands for not significant.

⁶Percent Pubmeth genes out of the prioritized candidate genes identified.

⁷Percent CGC genes out of the prioritized candidate genes identified.

We also compared the prioritized candidate genes that appear in more than one module selected. There are 68 prioritized candidate genes selected by more than one of the top 1% of modules identified by the proposed NEpiC algorithm, while there are 27 using the NEpi algorithm (Table 2 Part II). We then chose the top 68 genes with the strongest combined or mean signals using EpiC and Epi. According to Pubmeth, there are 8, 1, 2, 1 genes reported to be differently methylated in cancers out of the prioritized candidate genes identified using NEpiC, NEpi, EpiC and Epi, respectively. There are 14, 3, 2, 3 genes that have been causally implicated with cancers out of the prioritized candidate genes identified using NEpiC, NEpi, EpiC and Epi according to CGC, respectively (Table 2 Part II). The proposed NEpiC algorithm generates the highest percentages of reported differentially methylated genes in cancers and causally implicated cancer genes out of the prioritized candidate genes identified. Moreover, the percentages of differentially methylated genes in cancers according to Pubmeth and causally implicated cancer genes according to CGC out of the prioritized candidate genes are higher than that out of the candidate genes for both the NEpiC and NEpi algorithms, which agrees with the simulation results and suggests that the prioritization procedure using selection frequencies to further prioritize candidate genes might be a useful step for further selecting outcome-related genes. Between the list of eight genes reported to be differentially methylated in cancers according to Pubmeth and the list of 14 genes reported to be causally implicated in cancers according to CGC, there are two genes in common. Among the 12 causally implicated cancer genes identified by the CGC database only, five genes were also reported to be aberrantly methylated (two in breast cancers and three in other cancers, all were published after the Pubmeth database was generated), five genes were reported to be aberrantly expressed or mutated in breast cancer, and the remaining two genes were also reported to be aberrantly expressed or mutated in other cancers (Additional Table 1 for BRCA, Part II). We then conducted a gene set enrichment analysis of ‘pathway in cancer’ from KEGG among the prioritized candidate genes (Table 2 Part II). It shows that the proposed NEpiC algorithm achieves the most significant enrichment P-value while NEpi achieves the second most significant enrichment P-value. ‘Pathway in cancer’ was not enriched among candidate genes identified by EpiC and Epi algorithms.

We display the histogram of gene-level combined signal scores of the candidate genes in the top 1% of modules identified using NEpiC with the TCGA BRCA dataset in Supplementary Figure S3. About a quarter of the candidate genes identified have very high gene-level combined signal scores while there are also some candidate genes with rather small gene-level combined signal scores but which are connected to genes with high gene-level combined signal scores.

We further show in Figure 3 the original methylation measures of the CpG site with the largest site-level combined signal score (top row) and the site-level combined scores of all sites (bottom row) in three selected genes from the top first-ranked module identified by NEpiC with the TCGA BRCA data. We note that (i) tumor samples in general have bigger variation in methylation measures and this difference in variation is not driven by a few outliers; (ii) genes with large gene-level combined signal scores usually have a high proportion of sites with large site-level combined signal scores. Similar plots for the TCGA LIHC and the CUMC HCC data are shown in Supplementary Figures S1 and S2.

Figure 3.

Original methylation measures of the CpG site with the largest site-level combined signal score (top row) and site-level combined signal scores of all sites (bottom row) in three genes in the top first-ranked module identified with the TCGA BRCA data.

Among genes uniquely identified by NEpiC but not the other three methods in the TCGA BRCA dataset, we examined ZNF652 for illustration purpose. ZNF652 has a gene-level combined signal score of 9.02 and ranked #6,799 by gene-level combined signal scores. Thus, ZNF652 was not selected by EpiC. Moreover, ZNF652 is in the identified module using NEpiC with six other genes: CPLX2, NEUROG1, GFI1, PRDM14, ETS1 and RUNX1T1 (Figure 4), which ranked #7, #97, #251, #263, #360, and #665 by gene-level combined signal scores, respectively. CPLX2 and NEUROG1 were reported to be aberrantly methylated in cancers (30,31). On the other hand, using NEpi with mean signal scores, the highest ranked module with ZNF652 does not make to the top 1% of modules, thus was not selected by NEpi. This module has six genes: NEDD9, PIK3CA, ZBTB47, TCF3, CBFA2T3 and ZNF652, with the mean signal scores ranked #38, #167, #313, #1465, #2848 and #3651. Thus, ZNF652 was not selected by Epi either. Note that ZNF652 was previously reported to be associated with breast cancer (32).

Figure 4.

An example of identified module with ZNF652 uniquely identified by the NEpiC algorithm.

TCGA LIHC data and CUMC HCC data

After the same quality control steps, there are 229 700 CpG sites from 19,257 genes for 50 matched tumor and adjacent normal pairs in the TCGA LIHC data. Of those, 12 565 genes are also in the PPI network, which contains 115 964 edges. We then applied the proposed NEpiC algorithm using the PINA PPI network to the TCGA LIHC DNA methylation data.

Table 3 Part I displays the number of candidate genes identified in the TCGA LIHC data and the numbers (percentages) of reported differentially methylated genes in cancers according to Pubmeth (28) and causally implicated cancer genes according to CGC (26) out of the candidate genes identified. The proposed NEpiC algorithm clearly outperforms the other methods ignoring either biological network information or variance signals in DNA methylation. Between the list of 11 genes reported to be differentially methylated in cancers based on Pubmeth and the list of 18 genes reported to be causally implicated in cancers based on CGC, there are 6 genes in common. Among the 12 causally implicated cancer genes identified by the CGC database only, two genes were reported to be aberrantly methylated in cancers other than liver cancer (both were published after the Pubmed database was generated), seven genes were reported to be aberrantly expressed in liver cancer, and the remaining 3 genes were reported to be aberrantly expressed or mutated in other cancers. (Additional Table 2 for LIHC, Section I, Part I). In the pathway enrichment analysis, we selected eight core liver cancer pathways (p53 signaling pathway, cell cycle regulation pathway, TERT pathway, WNT pathway, chromatin modifying factors, growth factor signaling pathway, KEAP1–NFE2L2 pathway and NOTCH pathway) (25) and used the ‘pathway in cancer’ from KEGG as the benchmark to evaluate the performance of the four methods. The ‘pathway in cancer’, the ‘p53 signaling pathway’ and the ‘WNT signaling pathway’ are enriched among the candidate genes identified using the proposed NEpiC algorithm while only the ‘pathway in cancer’ is enriched among the candidate genes identified using NEpi (Table 3 Part I). All the top 1% of modules identified by NEpiC and NEpi with the TCGA LIHC data have permutation P-values smaller than 0.0005.

Number of candidate genes (Part I) and prioritized candidate genes (Part II) identified in the TCGA LIHC data and number of reported differentially methylated genes in cancers and causally implicated cancer genes out of the candidate genes (Part I) and out of the prioritized candidate genes (Part II) identified according to PubMeth¹ and CGC² and Bonferroni adjusted enrichment P-values of enriched core pathways in liver cancer and ‘pathway in cancer’ from KEGG among the candidate genes (Part I) and the prioritized candidate genes (Part II) identified

Table 3.

Number of candidate genes (Part I) and prioritized candidate genes (Part II) identified in the TCGA LIHC data and number of reported differentially methylated genes in cancers and causally implicated cancer genes out of the candidate genes (Part I) and out of the prioritized candidate genes (Part II) identified according to PubMeth¹ and CGC² and Bonferroni adjusted enrichment P-values of enriched core pathways in liver cancer and ‘pathway in cancer’ from KEGG among the candidate genes (Part I) and the prioritized candidate genes (Part II) identified

	NEpiC	NEpi	EpiC	Epi
Part I
Number of candidate genes identified	201	148	201	201
Number (percentage) of reported differentially methylated genes in cancers out of the candidate genes according to Pubmeth (%3)	11 (5.5%)	4 (2.7%)	8 (3.9%)	2 (0.1%)
Number (percentage) of causally implicated cancer out of the candidate genes according to CGC (%4)	18 (9.0%)	4 (2.7%)	10 (4.9%)	1 (0.5%)
Enrichment P-values5
Pathway in cancer	|$4.38 \times {10^{ - 10}}$|	|$1.78 \times {10^{ - 7}}$|	0.032	NS6
p53 signaling pathway	0.036	NS	NS	NS
WNT signaling pathway	|$6.49 \times {10^{ - 5}}$|	NS	NS	NS
Part II
Number of prioritized candidate genes identified	64	37	64	64
Number (percentage) of reported differentially methylated genes in cancers out of the prioritized candidate genes according to Pubmeth (%7)	5 (7.8%)	3 (8.1%)	3 (4.7%)	0 (0.0%)
Number (percentage) of causally implicated cancer out of the prioritized candidate genes according to CGC (%8)	10 (15.6%)	3 (8.1%)	4 (6.3%)	0 (0.0%)
Enrichment P-values9
Pathway in cancer	|$1.24 \times {10^{ - 8}}$|	|$1.80 \times {10^{ - 3}}$|	NS	NS
p53 signaling pathway	|$9.00 \times {10^{ - 4}}$|	NS	NS	NS
WNT signaling pathway	|$9.18 \times {10^{ - 7}}$|	NS	NS	NS

¹There are 292 reported differentially methylated genes in cancers according to Pubmeth (28).

²There are 572 mutated genes that have been causally implicated with cancers according to the Cancer Gene Census category (CGC, as of December 2015) (26).

³Percent Pubmeth genes out of the candidate genes identified.

⁴Percent CGC genes out of the candidate genes identified.

⁵Enrichment P-values of significant core liver cancer pathways and ‘pathway in cancer’ from KEGG among the candidate genes identified were Bonferroni corrected with the number of compared pathways from KEGG.

⁶NS stands for not significant.

⁷Percent Pubmeth genes out of the prioritized candidate genes identified.

⁸Percent CGC genes out of the prioritized candidate genes identified.

⁹Enrichment P-values of significant core liver cancer pathways and ‘pathway in cancer’ from KEGG among the prioritized candidate genes identified were Bonferroni corrected with the number of compared pathways from KEGG.

We then compared the prioritized candidate genes identified using the four methods (Table 3 Part II). There are 64 candidate genes selected by more than one of the top 1% of modules identified using NEpiC, while there are 37 using NEpi. We then chose the top 64 genes with the strongest combined or mean signals using EpiC and Epi algorithms. The percentages of reported differentially methylated genes in cancers according to Pubmeth and causally implicated cancer genes according to CGC among the prioritized candidate genes are higher than that among candidate genes before the prioritization procedure. This again suggests that the prioritization procedure may improve the performance of the proposed NEpiC algorithm. Between the list of five genes reported to be differentially methylated in cancers according to Pubmeth and the list of 10 genes reported to be causally implicated cancer genes according to CGC, there are three genes in common. Among the seven causally implicated cancer genes identified by the CGC database only, one gene was reported to be aberrantly methylated in cancer other than liver cancer (it was published after the Pubmeth database was generated), four genes were reported to be aberrantly expressed in liver cancer, and the remaining two genes were reported to be aberrantly expressed or mutated in other cancers (Additional Table 2 for LIHC, Section I, Part II). In the pathway enrichment analysis of the eight core liver cancer pathways and the ‘pathway in cancer’ from KEGG, three pathways (‘pathway in cancer’, ‘p53 signaling pathway’ and ‘WNT signaling pathway’) are enriched among the prioritized candidate genes identified by the proposed NEpiC algorithm, while only one pathway (‘pathway in cancer’) is enriched among the prioritized candidate genes identified by the NEpi algorithm (Table 3 Part II). No pathway is enriched among the prioritized candidate genes identified by NEpi and Epi algorithms.

We further performed a replication analysis using the CUMC HCC data and investigated the replication results using the TCGA LIHC data and the CUMC HCC data. Similar results as in Table 3 but using the CUMC HCC data are listed in Table S1 in the supplementary file. For replication analyses, we first define replicated candidate genes (replicated prioritized candidate genes) as the overlapping candidate genes (prioritized candidate genes) between the two lists of candidate genes (prioritized candidate genes) identified by the same method applied to the two HCC datasets. We then examined the reported differentially methylated genes in cancers according to Pubmeth and causally implicated cancer genes according to CGC out of the replicated candidate genes (replicated prioritized candidate genes) identified and the enrichment P-values of the replicated candidate genes (replicated prioritized candidate genes) in the eight liver cancer core pathways (25) and the ‘pathway in cancer’ from KEGG (Table 4). The proposed NEpiC algorithm generates the highest percentages of reported differentially methylated genes in cancers according to Pubmeth and causally implicated cancer genes according to CGC out of the replicated candidate genes and replicated candidate genes identified by the proposed NEpiC algorithm are enriched in two liver cancer core pathways and the ‘pathway in cancer’ (Table 4 Part I). Between the list of three genes reported to be differentially methylated in cancers according to Pubmeth and the list of five causally implicated genes according to CGC, there are two genes in common. Among the three causally implicated cancer genes identified by the CGC database only, one gene was also reported to be differentially methylated in cancer other than liver cancer (it was published after the Pubmeth database was generated), one gene was reported to be aberrantly expressed in liver cancer, and the last gene was reported to be aberrantly expressed in other cancer (Additional Table 2 for LIHC, Section II, Part I). For the replicated prioritized candidate genes, the proposed NEpiC algorithm also outperforms the other three methods using both percentages of causally implicated cancer genes according to CGC and enrichment P-values of eight liver cancer core pathways (25) and the ‘pathway in cancer’ from KEGG as the criteria (Table 4 Part II). However, the percentage of reported differentially methylated genes in cancers according to Pubmeth for NEpiC is lower than that for EpiC. This may due to the fact that the Pubmeth database is not most up-to-date. We also repeated the replication analyses with different cutoffs to select top ranked modules, where we chose candidate genes as the genes in the top 1.5%, 2%, 2.5% and 3% of modules (Supplementary Tables S2–S5). Under all these scenarios, NEpiC outperforms the other three methods which ignore either the biological network information or variance signals in both candidate genes and prioritized candidate genes using all three comparison criteria.

Number of replicated candidate genes (Part I) and replicated prioritized candidate genes (Part II) identified in the TCGA LIHC data and the CUMC HCC data and number of reported differentially methylated genes in cancers and causally implicated cancer genes out of the replicated candidate genes (Part I) and out of the replicated prioritized candidate genes (Part II) identified according to PubMeth¹ and CGC²and Bonferroni adjusted enrichment P-values of enriched core pathways in liver cancer and the ‘pathway in cancer’ from KEGG among the replicated candidate genes (Part I) and the replicated prioritized candidate genes (Part II) identified

Table 4.

Number of replicated candidate genes (Part I) and replicated prioritized candidate genes (Part II) identified in the TCGA LIHC data and the CUMC HCC data and number of reported differentially methylated genes in cancers and causally implicated cancer genes out of the replicated candidate genes (Part I) and out of the replicated prioritized candidate genes (Part II) identified according to PubMeth¹ and CGC²and Bonferroni adjusted enrichment P-values of enriched core pathways in liver cancer and the ‘pathway in cancer’ from KEGG among the replicated candidate genes (Part I) and the replicated prioritized candidate genes (Part II) identified

	NEpiC	NEpi	EpiC	Epi
Part I
Number of replicated candidate3 genes identified	42	26	89	77
Number (percentage) of reported differentially methylated genes in cancers out of the replicated candidate genes according to Pubmeth (%4)	3 (7.1%)	0 (0%)	3 (3.4%)	0 (0%)
Number (percentage) of causally implicated cancer out of the replicated candidate genes according to CGC (%5)	5 (11.9%)	0 (0%)	5 (5.6%)	0 (0%)
Enrichment P-values6
Pathway in cancer	|$2.56 \times {10^{ - 7}}$|	NS7	|$3.17 \times {10^{ - 4}}$|	NS
p53 signaling pathway	0.018	NS	NS	NS
WNT signaling pathway	|$3.60 \times {10^{ - 3}}$|	NS	0.028	NS
Part II
Number of replicated prioritized candidate genes identified	20	12	24	23
Number (percentage) of reported differentially methylated genes in cancers out of the replicated prioritized candidate genes according to Pubmeth (%8)	0 (0%)	0 (0%)	2 (8.3%)	0 (0%)
Number (percentage) of causally implicated cancer out of the replicated prioritized candidate genes according to CGC (%9)	1 (5.0%)	0 (0%)	1 (4.2%)	0 (0.0%)
Enrichment P-values10
p53 signaling pathway	0.031	NS	NS	NS

¹There are 292 reported differentially methylated genes in cancers according to Pubmeth (28).

²There are 572 mutated genes that have been causally implicated with cancers according to the Cancer Gene Census category (CGC, as of December 2015) (26).

³Replicated candidate genes are defined as the candidate genes identified in both the TCGA LIHC data and the CUMC HCC data.

⁴Percent Pubmeth genes out of the replicated candidate genes identified.

⁵Percent CGC genes out of the replicated candidate genes identified.

⁶Enrichment P-values of significant core liver cancer pathways and the ‘pathway in cancer’ from KEGG among the replicated candidate genes identified were Bonferroni corrected with the number of compared pathways from KEGG.

⁷NS stands for not significant.

⁸Percent Pubmeth genes out of the replicated prioritized candidate genes identified.

⁹Percent CGC genes out of the replicated prioritized candidate genes identified.

¹⁰Enrichment P-values of significant core liver cancer pathways and ‘pathway in cancer’ from KEGG among the replicated prioritized candidate genes identified were Bonferroni corrected with the number of compared pathways from KEGG.

A further investigation of the 42 replicated candidate genes from the top 1% of modules identified using the proposed NEpiC algorithm suggests that 33 replicated candidate genes out of the 42 replicated candidate genes have been reported to be associated with cancer. These includes genes are reported to be related with liver cancer: CDKN2A (33), GRASP (34), DLGAP1 (35), LPAR2 (35), STEAP4 (36), WNT3A (37), TSC22D1 (38), NKD2 (39), TGFA (40), TERT (41), HSP90AA1 (42) and IGF1R (43); genes that are known to be aberrantly methylated in cancers other than liver cancer: MYO10 (44), BAI1 (45), FYN (46), ACTA1 (47), SPRR2A (48) and CARD11 (49); and genes that are associated with cancers other than liver cancer: VIM (50), CTBP2 (51), PRKCQ (52), PDZD2 (53), DSCAML1 (54), KCNQ1 (55), KCNQ2 (56), KCNQ3 (57), OBSCN (58), FSCN1 (59), DLGAP2 (60), CFTR (61), RUNX1T1 (62), GRID2 (63) and SCN5A (64). The full list of 42 replicated candidate genes is included in the supplementary materials.

DISCUSSION

In this article, we proposed the NEpiC algorithm, a network assisted method incorporating combined signals in mean and variance differences of DNA methylation data. We demonstrated that incorporating prior biological network information and utilizing the signals in variance differences of DNA methylation data could effectively improve the power of the association studies to identify aberrant methylated genes associated with the outcomes.

We demonstrated a much improved power of the proposed NEpiC algorithm that incorporates both biological network information and variance signals in DNA methylation data than the methods that do not. In simulation studies, the proposed NEpiC algorithm identifies most truly associated genes among the candidate genes identified and achieves the most significant enrichment P-value of truly associated genes among candidate genes identified. Using the prioritized genes that were selected in more than one of the top 1% of modules further improves the performance of the proposed NEpiC algorithm in both simulation studies and real data applications. The application to two independent liver cancer datasets, the TCGA LIHC data and the CUMC HCC data, gives us the opportunity to examine replication results. The replication results show that the proposed NEpiC algorithm identifies more genes that were already reported to be differentially methylated in cancers according to Pubmeth and causally implicated cancer genes according to CGC and identifies genes that are more enriched in known liver cancer pathways than methods that do not use both biological network information and variance signals in DNA methylation data.

Although we focused on applying the proposed NEpiC algorithm to cancer patients with tumor and adjacent normal tissues to identify genes that are related to tumor status, several publications (65,66) have shown that differential variability is most informative and meaningful in comparing precursor cancer lesions to normal cells. That is, the developed NEpiC algorithm may be the most useful in cancer early detection.

Since bigger variation in methylation measures is usually observed in tumor tissues compared to adjacent tissues, tumor tissue purity may thus influence findings from methods that use variance signals. A method was recently developed to check the purity of tumor tissues using DNA methylation 450K arrays (67), which could be applied in the quality control steps to screen out tumor tissues with low purity.

In summary, we developed a new algorithm, NEpiC that incorporates biological network information and utilizes variance signals in DNA methylation data in detecting differentially methylated genes. Results from both simulations and real data applications demonstrate a much better performance of the NEpiC algorithm compared to several other methods that ignore either the biological network information or variance signals in DNA methylation data. The NEpiC algorithm is implemented in an R package which is freely available through CRAN.

FUNDING

National Natural Science Foundation of China [61272380 to P.R. and S. Z.]; National Institutes of Health [R01ES05116 and P30ES009089 to S.J., R.S. and S.W.]. Funding for open access charge: Departmental fund from Department of Biostatistics, Columbia University.

Conflict of interest statement. None declared.