DSLE2 random-effects meta-analysis model for high-throughput methylation data

Abstract

Background

With the rapid development of high-throughput sequencing technology, high-throughput sequencing data has grown on a massive scale, leading to the emergence of multiple public databases, such as EBI and GEO. Conducting secondary mining of high-throughput sequencing data in these databases can yield more valuable insights. Meta-analysis can quantitatively combine high-throughput sequencing data from the the same topic. It increases the sample size for data analysis, enhances statistical power, and results in more consistent and reliable conclusions.

Results

This study proposes a new between-study variance estimator $E_{m}$. We prove that $E_{m}$ is non-negative and $E_{m} \left( {\hat{\tau }_{m}^{2} } \right)$ increases with the increase of $\hat{\tau }_{m}^{2}$, satisfying the general conditions of the between-study variance estimator. We get the DSLE2 (two-step estimation starting with the DSL estimate and the $E_{m}$ in the second step) random-effects meta-analysis model based on the between-study variance estimator Em. The accuracy and a series of evaluation metrics of the DSLE2 model are better than those of the other 6 meta-analysis models. DSLE2 model is applied to lung cancer and Parkinson’s methylation data. Significantly differentially methylated sites identified by DSLE2 model and the genes with significantly differentially methylated sites are closely related to two diseases, indicating the effectiveness of DSLE2 random-effects model.

Conclusions

This paper propose the DSLE2 random-effects meta-analysis model based on new between-study variance estimator Em. The DSLE2 model performs well for methylation data.

Background

DNA methylation refers to the addition of a methyl group at the CpG site of the DNA chain [1]. DNA methylation occurs at millions of CpG dinucleotide positions in the genome and changes with age and the external environment [2]. It is an important epigenetic modification [3]. In the correlation analysis between methylation and traits, HM450K (Infinium Human Methylation 450 K BeadChip) or other Illumina technologies use the methylated fluorescence signal intensity ($M$) and the unmethylated fluorescence signal intensity ($U$) of the CpG site to calculate DNA methylation levels [4]. Two commonly used calculation methods for methylation levels are $\beta = M/(M + U)$ and $M = \log_{2} (M/U)$ [5]. In the minfi R package, the getBeta function and getM function can be used directly to calculate the $\beta$ value and $M$ value of every probe [6].

With the development of high-throughput sequencing technology, multiple public databases have been formed, such as the GEO (Gene Expression Omnibus) database and TCGA (The Cancer Genome Atlas database) database [7]. People can obtain methylation sequencing data of the same disease from different databases [8]. Due to differences in experimental conditions and sample processing procedures, the results of methylation-trait association analysis of the same disease in different studies are different [9]. Meta-analysis can combine methylation data on the same topic to provide a more reliable list of significantly different methylation sites [10]. In meta-analysis, datasets from independent studies with the same purpose often have heterogeneity [11]. In the random-effects model, the comprehensive effect-size result depends on the quantification of heterogeneity between studies [12]. Therefore, the estimation of heterogeneity is the most critical step in the random-effects meta-analysis model [13, 14]. In addition to quantifying heterogeneity, investigating differences in experimental designs, participants, and interventions across studies can help to understand the causes of heterogeneity [15, 16]. This paper aims to estimate the size of the between-study variance in a random-effects meta-analysis model.

We propose the DSLE2 meta-analysis model based on new between-study variance estimator $E_{m}$ for methylation high-throughput sequencing data. Under three hypothesis testing conditions, this model is compared with DSLR2 (two-step estimation starting with the DSL estimate and the $R^{2}$ in the second step), DL (DerSimonian and Laird estimate), EB (Empirical Bayes estimate), HO (Hedges and Olkin estimate), RML (restricted maximum likelihood estimate), and SJ (Sidik and Jonkman estimate) random-effects meta-analysis models using sensitivity, and other evaluation metrics. The results show that the DSLE2 methylation meta-analysis model performs well under the first hypothesis testing condition. We apply the DSLE2 random-effects model to lung cancer and Parkinson’s methylation data, further demonstrating the reliability of the DSLE2 meta-analysis model.

The proposed between-study variance estimator is designed to accurately quantify the true variability in effect sizes across studies, distinguishing it from random error. This advancement enhances the precision and reliability of meta-analysis models, offering significant benefits for research methodology, personalized medicine, study design, and resource allocation. By implementing and validating this new estimator, we can achieve a clearer understanding of effect consistency across studies, ultimately leading to more robust and generalizable conclusions. The primary purpose of the proposed between-study variance estimator is to improve the modeling of heterogeneity in methylation high-throughput sequencing data. FEM meta-analysis models often struggle to accurately capture the complex variability inherent in such data due to differences in study design, populations, and technical procedures. Accurate effect size estimation is critical for identifying biologically significant methylation patterns and their potential clinical implications. To refine the estimation of true effect sizes by more accurately accounting for between-study variance is very important. Reliable identification of consistent methylation markers across studies is also crucial for developing personalized medical interventions. The improved model supports the discovery of biomarkers that are consistently reproducible, aiding in disease diagnosis, prognosis, and treatment customization.

Methods

The random-effects model is a generalization of the fixed-effects model. The fixed effects model assumes that all studies included in the same meta-analysis have the same true effect size, while the random-effects model assumes that the true effect sizes of different studies included in the same meta-analysis obey a normal distribution. The between-study variance is a statistic used by the random-effects model to measure the heterogeneity between studies from the same topic. If the between study variance is 0, the random-effects model degenerates into the fixed-effects model. The between-study variance estimator of DSL (DerSimonian and Laird estimate) random-effects model is simple and the most commonly used method for estimating the between-study heterogeneity. In addition to being easy to calculate, the between-study variance estimator is also suitable for effect sizes of different dimensions. However, in practice, due to the possible occurrence of negative values, DSL between-study variance estimators are often truncated. In this paper, we present a non-truncated estimator of between-study variance $E_{m}$.

DSLE2 random-effects meta-analysis model

We propose a general heterogeneity variance estimator for methylation sequencing data that is applicable to effect sizes at any scale and is non-negative. Assume that $y_{1m} ,y_{2m} , \cdots ,y_{km}$ are the effect sizes of k independent studies of methylation site m; $Q_{m}$ is a heterogeneity statistic for methylation site m; $\omega_{im}$ is the weight of methylation site m for ith study in the fixed-effects meta-analysis model; $\omega_{im}^{*}$ is the weight of methylation site m for ith study in the random-effects meta-analysis model; $\sigma_{im}^{2}$ is the within-study variability representing sampling errors of ith study; $\hat{\tau }_{m,DSL}^{2}$ is the between-study variance estimator in DSL random-effects model; $\hat{\mu }_{m,DSL}$ is the mean of effects in DSL random-effects model. And the between-study variance estimator $E_{m}$ is

$$E_{m} = \frac{{Q_{m} - S_{MM,m} }}{{\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } - \sum\limits_{i = 1}^{k} {(\omega_{im}^{*} )^{2} } /\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } }}$$

(1)

where

$$Q_{m} = \sum\limits_{i = 1}^{k} {\omega_{im} } \left( {y_{im} - \hat{\mu }_{m} } \right)^{2}$$

$$\hat{\mu }_{m} = \frac{{\sum\limits_{i = 1}^{k} {\omega_{im} } y_{im} }}{{\sum\limits_{i = 1}^{k} {\omega_{im} } }}$$

$$\omega_{im} = \frac{1}{{\sigma_{im}^{2} }},i = 1,2, \cdots ,k$$

$$S_{MM,m} = \sum\limits_{i = 1}^{k} {\omega_{im}^{*} } \left( {y_{im} - \hat{\mu }_{{m,{\text{DSL}}}} } \right)^{2}$$

$$\hat{\mu }_{{m,{\text{DSL}}}} = \frac{{\sum\limits_{i = 2}^{k} {\omega_{im}^{*} } y_{im} }}{{\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } }}$$

$$\omega_{im}^{*} = \frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{{m,{\text{DSL}}}}^{2} }},i = 1,2, \cdots ,k$$

The algorithm of DLSE2 meta-analysis model is as follows:

First, we use the fixed-effects model to calculate the weight and comprehensive effect value of each study

$$\omega_{im} = \frac{1}{{\sigma_{im}^{2} }},i = 1,2, \cdots ,k$$

and

$$\hat{\mu }_{m} = \frac{{\sum\limits_{i = 1}^{k} {\omega_{im} } y_{im} }}{{\sum\limits_{i = 1}^{k} {\omega_{im} } }}$$

Then, we calculate the heterogeneity statistic $Q_{m}$

$$Q_{m} = \sum\limits_{i = 1}^{k} {\omega_{im} } (y_{im} - \hat{\mu }_{m} )^{2}$$

We use the DSL random-effects model to calculate the between-study variance estimator, the weight of each study, and the corresponding comprehensive effect size

$$\hat{\tau }_{{m,{\text{DSL}}}}^{2} = \max \left\{ {0,\frac{{Q_{m} - (k - 1)}}{{S_{1m} - S_{2m} /S_{1m} }}} \right\}$$

$$S_{rm} = \sum\limits_{i = 1}^{k} {\omega_{im}^{r} } ,r = 1,2$$

$$\omega_{im}^{*} = \frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{{m,{\text{DSL}}}}^{2} }},i = 1,2, \cdots ,k$$

$$\hat{\mu }_{{m,{\text{DSL}}}} = \frac{{\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } y_{im} }}{{\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } }}$$

We further calculate $S_{MM,m}$

$$S_{MM,m} = \sum\limits_{i = 1}^{k} {\omega_{im}^{*} } \left( {y_{im} - \hat{\mu }_{{m,{\text{DSL}}}} } \right)^{2}$$

then

$$E_{m} = \frac{{Q_{m} - S_{MM,m} }}{{\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } - \sum\limits_{i = 1}^{k} {\omega_{im}^{*2} } /\sum\limits_{i = 1}^{k} {\omega_{im}^{*} } }}$$

The weight of each study based on the between-study variance $E_{m}$ is

$$\omega_{{im,{\text{DSLE2}}}} = \frac{1}{{\sigma_{im}^{2} + E_{m} }},i = 1,2, \cdots ,k$$

The comprehensive effect is

$$\hat{\mu }_{{m,{\text{DSLE2}}}} = \frac{{\sum\limits_{i = 1}^{k} {\omega_{{im,{\text{DSLE2}}}} } y_{im} }}{{\sum\limits_{i = 1}^{k} {\omega_{{im,{\text{DSLE2}}}} } }}$$

The random-effect for each study is

$$\hat{\xi }_{ig} = \frac{{E_{m} }}{{\sigma_{im}^{2} + E_{m} }}\left( {y_{im} - \hat{\mu }_{{m,{\text{DSLE2}}}} (\hat{\tau }_{{m,{\text{DSLE2}}}}^{2} )} \right),i = 1,2, \cdots ,k$$

The variance of comprehensive effect is

$${\text{Var}}(\hat{\mu }_{{m,{\text{DSLE2}}}} ) = \frac{1}{{\sum\limits_{i = 1}^{k} {\omega_{{im,{\text{DSLE2}}}} } }}$$

The $z$ statistic of comprehensive effect is

$$z_{m} (\hat{\mu }_{{m,{\text{DSLE2}}}} ) = \frac{{\hat{\mu }_{{m,{\text{DSLE2}}}} }}{{\sqrt {{\text{Var}}(\hat{\mu }_{{m,{\text{DSLE2}}}} )} }}\sim N(0,1)$$

The lower and upper bounds of the $100(1 - \alpha )\%$ confidence interval are

$$LL\left[ {\hat{\mu }_{{m,{\text{DSLE2}}}} } \right] = \hat{\mu }_{{m,{\text{DSLE2}}}} - z_{(1 - \alpha )/2(m)} \times \sqrt {{\text{Var}}\left[ {\hat{\mu }_{{m,{\text{DSLE2}}}} } \right]}$$

and

$$UL\left[ {\hat{\mu }_{{m,{\text{DSLE2}}}} } \right] = \hat{\mu }_{{m,{\text{DSLE2}}}} + z_{(1 - \alpha )/2(m)} \times \sqrt {{\text{Var}}\left[ {\hat{\mu }_{{m,{\text{DSLE2}}}} } \right]}$$

the $p$-value of one-sided test is

$$p = 1 - \Phi \left( {z_{m} \left( {\hat{\mu }_{{m,{\text{DSLE2}}}} } \right)} \right)$$

the $p$-value of two-sided test is

$$p = 2\left[ {1 - \Phi \left( {|z_{m} \left( {\hat{\mu }_{{m,{\text{DSLE2}}}} |} \right)} \right)} \right]$$

where $\Phi ( \cdot )$ is the cumulative function of the standard normal distribution.

Theorem 1 Assume that $y_{1m} ,y_{2m} , \cdots ,y_{nm}$ are the effect sizes of n independent studies of methylation site $m$, and the between-study variance estimator $E_{m}$ is

$$E_{m} = \frac{{Q_{m} - S_{MM,m} }}{A}$$

(2)

where

$$Q_{m} = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} \left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} y_{i} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} }}} \right)^{2}$$

(3)

$$S_{MM,m} = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} y_{im} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}} \right)^{2} \,$$

(4)

$$A = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} - \frac{{\sum\limits_{i = 1}^{n} {\left( {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \right)^{2} } }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}$$

(5)

Then,

$E_{m} (\hat{\tau }_{m}^{2} )$ is monotone and non-decreasing with respect to $\hat{\tau }_{m}^{2}$.

Proof It can be obtained from (3)

$$\frac{{dQ_{m} }}{{d\hat{\tau }_{m}^{2} }} = 0$$

(6)

It can be obtained from (4)

$$\frac{{dS_{MM,m} (\hat{\tau }_{m}^{2} )}}{{d\hat{\tau }_{m}^{2} }} = - \sum\limits_{i = 1}^{n} {\left( {\left( {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \right)\left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} y_{im} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}} \right)} \right)^{2} } \le 0$$

(7)

It can be obtained from (5)

$$\frac{{dA(\hat{\tau }_{m}^{2} )}}{{d\widehat{{\tau_{m}^{2} }}}} = - \sum\limits_{i = 1}^{n} {\frac{1}{{(\sigma_{im}^{2} + \hat{\tau }_{m}^{2} )^{2} }}} + \frac{{2\sum\limits_{i = 1}^{n} {\frac{1}{{(\sigma_{im}^{2} + \hat{\tau }_{m}^{2} )^{3} }}} \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} - \left( {\sum\limits_{i = 1}^{n} {\frac{1}{{(\sigma_{im}^{2} + \hat{\tau }_{m}^{2} )^{2} }}} } \right)^{2} }}{{\left( {\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} } \right)^{2} }}$$

let

$$x_{i} = \frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}$$

Then

$$x_{i}> 0$$

$$\frac{{{\text{d}}A(\hat{\tau }_{m}^{2} )}}{{{\text{d}}\hat{\tau }_{m}^{2} }} = - \sum\limits_{i = 1}^{n} {x_{i}^{2} } + \frac{{2\sum\limits_{i = 1}^{n} {x_{i}^{3} } \sum\limits_{i = 1}^{n} {x_{i} } - (\sum\limits_{i = 1}^{n} {x_{i}^{2} } )^{2} }}{{\left( {\sum\limits_{i = 1}^{n} {x_{i} } } \right)^{2} }}$$

$$= \frac{{ - \sum\limits_{i = 1}^{n} {x_{i}^{2} } \left( {\sum\limits_{i = 1}^{n} {x_{i} } } \right)^{2} + 2\sum\limits_{i = 1}^{n} {x_{i}^{3} } \sum\limits_{i = 1}^{n} {x_{i} } - \left( {\sum\limits_{i = 1}^{n} {x_{i}^{2} } } \right)^{2} }}{{\left( {\sum\limits_{i = 1}^{n} {x_{i} } } \right)^{2} }}$$

$$\le \frac{{2\left( {\sum\limits_{i = 1}^{n} {x_{i}^{3} } \sum\limits_{i = 1}^{n} {x_{i} } - \left( {\sum\limits_{i = 1}^{n} {x_{i}^{2} } } \right)^{2} } \right)}}{{\left( {\sum\limits_{i = 1}^{n} {x_{i} } } \right)^{2} }}$$

Let

$$f(y) = \sum\limits_{i = 1}^{n} {x_{i}^{(4 - y)} } \sum\limits_{i = 1}^{n} {x_{i}^{y} }$$

Then

$$\frac{{{\text{d}}f}}{{{\text{d}}y}} = - \sum\limits_{i = 1}^{n} {x_{i}^{(4 - y)} } \ln x_{i} \sum\limits_{i = 1}^{n} {x_{i}^{y} } + \sum\limits_{i = 1}^{n} {x_{i}^{(4 - y)} } \sum\limits_{i = 1}^{n} {x_{i}^{y} } \ln x_{i}$$

So, when $y = 2$, $\frac{{{\text{d}}f}}{{{\text{d}}y}} = 0$.

We can get that $y = 2$ is the stagnation point of $f(y)$.

Then, we can get

$$\sum\limits_{i = 1}^{n} {x_{i}^{3} } \sum\limits_{i = 1}^{n} {x_{i} } \le \left( {\sum\limits_{i = 1}^{n} {x_{i}^{2} } } \right)^{2}$$

So

$$\frac{{{\text{d}}A(\hat{\tau }_{m}^{2} )}}{{{\text{d}}\hat{\tau }_{m}^{2} }} \, \le \frac{{2\left( {\sum\limits_{i = 1}^{n} {x_{i}^{3} } \sum\limits_{i = 1}^{n} {x_{i} } - \left( {\sum\limits_{i = 1}^{n} {x_{i}^{2} } } \right)^{2} } \right)}}{{\left( {\sum\limits_{i = 1}^{n} {x_{i} } } \right)^{2} }} \le 0$$

(8)

It can be obtained from (2)

$$\frac{{{\text{d}}(E_{m} (\hat{\tau }_{m}^{2} ))}}{{{\text{d}}\hat{\tau }_{m}^{2} }} = \frac{{{\text{d}}\left( {\frac{{Q_{m} - S_{MM,m} (\hat{\tau }_{m}^{2} )}}{{A(\hat{\tau }_{m}^{2} )}}} \right)}}{{{\text{d}}\hat{\tau }_{m}^{2} }} = \frac{{ - \left( {A(\hat{\tau }_{m}^{2} )\frac{{{\text{d}}S_{MM,m} (\hat{\tau }_{m}^{2} )}}{{{\text{d}}\hat{\tau }_{m}^{2} }} + \left( {Q_{m} - S_{MM,m} } \right)\frac{{{\text{d}}A(\hat{\tau }_{m}^{2} )}}{{{\text{d}}\hat{\tau }_{m}^{2} }}} \right)}}{{\left( {A(\hat{\tau }_{m}^{2} )} \right)^{2} }}$$

(9)

It can be obtained from (3) and (4)

$$Q_{m} = S_{MM,m} (0) \,$$

It is attainable via (7)

$$Q_{m} - S_{MM,m} (\hat{\tau }_{m}^{2} ) \ge 0$$

(10)

It is available through (6), (7), (8) and (10)

$$\frac{{{\text{d}}(E_{m} (\hat{\tau }_{m}^{2} ))}}{{{\text{d}}\hat{\tau }_{m}^{2} }} \ge 0 \,$$

So, $E_{m} (\hat{\tau }_{m}^{2} )$ increases with the increase of $\hat{\tau }_{m}^{2}$.

Theorem 2 Assume that $y_{1m} ,y_{2m} , \cdots ,y_{nm}$ are the effect sizes of n independent studies of methylation site m, and the between-study variance estimator ${\text{E}}_{{\text{m}}}$ is

$$E_{m} = \frac{{Q_{m} - S_{MM,m} }}{A}$$

where

$$Q_{m} = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} \left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} y_{i} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} }}} \right)^{2}$$

$$S_{MM,m} = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} y_{im} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}} \right)^{2}$$

$$A = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} - \frac{{\sum\limits_{i = 1}^{n} {\left( {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \right)^{2} } }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}$$

Then, $E_{m} \ge 0$.

Proof Because $A(\hat{\tau }_{m}^{2} ) \, = \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} - \frac{{\sum\limits_{i = 1}^{n} {\left( {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \right)^{2} } }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}> 0$

$$Q_{m} - S_{MM,m} (\hat{\tau }_{m}^{2} ) \ge 0$$

we have

$$E_{m} (\hat{\tau }_{m}^{2} ) \, = \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} \left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} y_{i} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} }}} }}} \right)^{2} - \sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \left( {y_{im} - \frac{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} y_{im} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}} \right)^{2} }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} - \frac{{\sum\limits_{i = 1}^{n} {\left( {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} \right)^{2} } }}{{\sum\limits_{i = 1}^{n} {\frac{1}{{\sigma_{im}^{2} + \hat{\tau }_{m}^{2} }}} }}}} = \frac{{Q_{m} - S_{MM,m} (\hat{\tau }_{m}^{2} )}}{{A(\hat{\tau }_{m}^{2} )}} \ge 0$$

Results

Simulation study

We simulated five sets of meta-analysis data. Each set of data has $K(K = 4,6,8,10,12)$ studies. Each study contains 20,000 methylation sites and 30 sample data. The first 15 samples are control data, and the last 15 sample data are experimental data. The methylation sites of each sample in each study are divided into 100 categories $(c = 1,2, \cdots ,100)$, with 200 methylation sites in each category. The top 9000 methylation sites of each study are divided into K groups $(k_{m} = 1,2, \cdots ,K)$. The first $9000/K$ methylation sites belong to the first group $(k_{m} = 1)$, and the $9000/K + 1$ to $9000*2/K$ methylation sites belong to the second group group $(k_{m} = 2)$, and so on, the methylation site from $9000(K - 1)/K + 1$ to 9000 methylation sites belong to the K-th group $(k_{g} = K)$, the 9001 to 20,000 methylation sites belong to the 0th group $(k_{m} = 0)$. The algorithm for simulating data is summarized as follows:

First, we randomly sampled the methylation sites of the $c(1 \le c \le 100)$ th class of the $k(1 \le k \le K)$ th study $k(1 \le k \le K)$, where $k(1 \le k \le K)$, $I_{200 \times 200}$ is identity matrix, $J_{200 \times 200}$ is a matrix with all elements $1$, $W^{ - 1}$ is the inverse of the Wishart distribution, and $\sum^{\prime}_{ck}$ is normalized to $\sum_{ck}$, such that all its diagonal elements are $1$. Then, sample the methylation level of the $c(1 \le c \le 100)$ th class site of the $n(1 \le n \le 30)$ sample of the $k$ th study $(X^{\prime}_{mc1nk} ,X^{\prime}_{mc2nk} , \cdots ,X^{\prime}_{mc200nk} )^{T} \sim MVN(0,\sum\limits{_{ck} )} ,1 \le k \le K$.

We performed differential methylation settings on the first 9000 methylation sites, randomly sampling $\delta_{mk} \in (0,1)$, so that $\sum\limits_{k = 1}^{K} {\delta_{mk} } = k_{m} (k_{m} = 1,2, \cdots ,K)$. When $\delta_{mk} = 1$, site $m$ of the $k$ th study is a significantly differentially methylated site. When $\delta_{mk} = 0$, site $m$ of the $k$ th study is a non-significantly differential methylation site. This paper conducts random sampling $\mu_{mk} \sim U(0.5,3)$. The methylation level of the control group remains unchanged, and the methylation expression level of the experimental group is $Y_{mnk} = X^{\prime}_{m(n + N)k} + \mu_{mk} *\delta_{mk}$, where $1 \le m \le 9000,1 \le n \le 15,1 \le k \le K$.

We compared the distribution of between study variance estimators of the DSLE2 random-effects model and six other meta-analysis models. As can be seen from Fig. 1, compared with the other six random-effects models, the between-study variance estimators of the DSLE2 meta-analysis model is relatively concentrated. The between-study variance estimators of DSLE2, DL, EB, HO, RML, and SJ random-effects models are mostly below 0.5. Compared with the other 6 random-effects models, the DSLE2 random-effects model has a relatively larger number of between-study variances distributed around 0.

Fig. 1

Violin plot of between-study variance for seven random-effects meta-analysis models

We also compared the accuracy, false negative rate, negative prediction rate, recall rate, Matthews correlation coefficient, PCMiss (Prediction-conditioned miss) value, SAR value, and F value of the DSLE2 meta-analysis model with DL, EB, HO, RML, and SJ meta-analysis models under three hypothesis.

Accuracy

We compared the accuracy of seven random-effects models under three hypothesis (see Fig. 2, Figure S1 and Figure S2). Under the first hypothesis, the DSLE2 meta-analysis model has the highest accuracy, the SJ random-effects model has the lowest accuracy, and the DSLR2, DL, EB, HO, and RML meta-analysis models have the similar accuracy. The accuracy of the DSLR2, DL, EB, HO, RML, and SJ meta-analysis models decreases as the number of studies increases. There is no obvious decreasing trend as the number increases for DSLE2 random-effects model. Under the second hypothesis, the DSLE2 meta-analysis model has the lowest accuracy, the SJ meta-analysis model has the highest accuracy.The accuracy of the DSLR2, DL, EB, HO, RML, and SJ meta-analysis models increases with the number of studies increases. Under the third hypothesis, the DSLE2 random-effects model has the lowest accuracy, and the DSLR2 meta-analysis model has the highest accuracy.

Fig. 2

Plot of accuracy for 7 random-effects meta-analysis models under the first hypothesis

False negative rate

The false negative rate, also known as the second type error rate, refers to the proportion of the number of significantly differential methylation sites predicted by the model as non-significantly differential methylation sites to the number of all significantly differential methylation sites. We compared the false negative rate of seven meta-analysis models under three hypothesis (see Fig. 3, Figure S3 and Figure S4).

Fig. 3

Plot of false negative rate curves for seven random-effects meta-analysis models under the first hypothesis

Under the first hypothesis, the false negative rate of the DSLE2 model is the lowest among the seven random-effects models; the false negative rate of the SJ model is the highest among the seven random-effects models. The false negative rates of the DSLR2, DL, EB, HO, and RML meta-analysis models are relatively close, and they are all lower than the false-negative rates of the SJ meta-analysis model, and they are all higher than the false negative rate of the DSLE2 meta-analysis model. The false negative rate of the DSLR2, DL, EB, HO, RML, and SJ meta-analysis models increases with the increase of the number of studies. Under the second hypothesis, the false negative rate of the DSLE2 meta-analysis model is the lowest among the seven random-effects models, and the false negative rates of the DL, EB, HO, and RML random-effects models are relatively close, and the false negative rate of the SJ model is the highest among the seven meta-analysis models.

Under the third hypothesis, the false negative rate of the DSLE2 meta-analysis model is the lowest among the seven meta-analysis models; the false negative rate of the SJ random-effects model is the highest among the seven random-effects models. The false negative rates of the DSLR2, DL, EB, HO, and RML meta-analysis models are relatively close, and they are all lower than the false negative rates of the SJ random-effects model, and they are all higher than the false negative rate of the DSLE2 meta-analysis model. The false negative rate of the DSLR2, DL, EB, HO, RML, and SJ random-effects models increases with the increase of the number of studies.

Matthews correlation coefficient

Matthews correlation coefficient is the correlation coefficient that describes the actual class and the predicted class. We compared the Matthews correlation coefficients of seven random-effects models under three hypothesis testing conditions (see Fig. 4, Figure S5 and Figure S6). Under the first hypothesis, when the number of studies is greater than 4, the Matthews correlation coefficient of the DSLE2 meta-analysis model is the highest among the Matthews correlation coefficient of the seven random-effects models. The Matthews correlation coefficient of SJ meta-analysis model is the lowest among the Matthews correlation coefficients of the seven meta-analysis models. The Matthews correlation coefficients of DSLR2, DL, EB, HO, and RML random-effects models are close, which are lower than the Matthews correlation coefficient of the DSLE2 meta-analysis model and higher than the Matthews correlation coefficient of the SJ random-effects model. The Matthews correlation coefficient of DSLR2, DL, EB, HO, RML, SJ meta-analysis models decreases as the number of studies increases. Under the second hypothesis, the Matthews correlation coefficient of the DSLE2 random-effects model is the lowest among the Matthews correlation coefficients of the seven meta-analysis models; the Matthews correlation coefficient of the SJ random-effects model is the highest among the Matthews correlation coefficients of the seven random-effects models; the Matthews correlation coefficients of DSLR2, DL, EB, HO, and RML meta-analysis models are relatively close. The Matthews correlation coefficients of HO, RML, and SJ random-effects models increase with the increase of the number of studies. Under the third hypothesis, the Matthews correlation coefficient of the DSLE2 random effects meta-analysis model is the lowest among the Matthews correlation coefficients of the seven meta-analysis models, and the Matthews correlation coefficient of the DSLR2 meta-analysis model is the highest among the Matthews correlation coefficients of the seven random-effects meta-analysis models.

Fig. 4

Plot of Matthews correlation coefficient curves for 7 random-effects models under the first hypothesis

Negative predictive value

Negative predictive value refers to the proportion of the number of non-significantly differentially methylated sites correctly predicted by meta-analysis model to the number of all non-significantly differentially methylated sites predicted by meta-analysis model. We compared the negative predictive values of seven meta-analysis models (see Fig. 5, Figure S7 and Figure S8). Under the first hypothesis, the negative predictive values of the DSLE2 meta-analysis model are the highest among seven random-effects models; the negative predictive values of the SJ meta-analysis model are the lowest among seven random-effects models. The negative predictive values of the DSLR2, DL, EB, HO, and RML random-effects models are relatively close, and they are all lower than the negative predictive value of the DSLE2 random-effects model and higher than the negative predictive values of SJ random-effects model. Under the second hypothesis, the negative predictive values of the SJ meta-analysis model are the lowest among the seven random-effects models. The negative predictive values of the DL, EB, HO, and RML random-effects models are relatively close. Under the third hypothesis, the negative predictive values of the DSLE2 meta-analysis model are the highest among the seven random effects models; the negative predictive values of the SJ random-effects model are the lowest among the seven meta-analysis models. The negative predictive rates of DSLR2, DL, EB, HO, and RML meta-analysis models are relatively close.

Fig. 5

Plot of negative predictive value curves for 7 random-effects models under the first hypothesis

Prediction-conditioned miss

PCMiss (Prediction-Conditioned Miss) refers to the proportion of the number of sites incorrectly predicted as non-significantly differentially methylated to the number of all non-significantly differentially methylated sites predicted by the meta-analysis model. We compared the PCMiss values of seven meta-analysis models under three hypothesis (see Fig. 6, Figure S9 and Figure S10). Under the first hypothesis, the PCMiss values of the DSLE2 random-effects model are the lowest among the PCMiss values of the seven meta-analysis models; the PCMiss values of the SJ random-effects model are the highest among the seven meta-analysis models. The PCMiss values of the DSLR2, DL, EB, HO, and RML random-effects meta-analysis models are relatively close; the PCMiss values of the seven models all increase with the increase of the number of studies. Under the second hypothesis, the PCMiss values of the SJ model are the highest among the PCMiss values of the seven models. The PCMiss values of the DL, EB, HO, and RML random-effects models are close. Under the third hypothesis testing, the PCMiss values of the DSLE2 model are the lowest among the PCMiss values of seven random effects meta-analysis models; the PCMiss values of the SJ meta-analysis model are the highest among the seven random-effects models. The PCMiss values of the DSLR2, DL, EB, HO, and RML meta-analysis models are relatively close.

Fig. 6

Plot of PCMiss-value curves for 7 random-effects models under the first hypothesis

Recall

This paper measured the recall rates of seven random-effects models under three hypothesis (see Fig. 7, Figure S11 and Figure S12). Under the first hypothesis, the recall rates of the DSLE2 meta-analysis model are the highest among seven random effects meta-analysis models; the recall rates of the SJ random-effects model are the lowest among seven meta-analysis models. The recall rates of the DL, EB, HO, and RML meta-analysis models are relatively close. The recall rates of the DSLR2, DL, EB, HO, RML, and SJ random-effects models decreases as the number of studies increases. Under the second hypothesis, the recall rates of the DSLE2 model are the highest among the seven meta-analysis models, followed by the recall rates of the DSLR2 model, and the recall rates of the SJ are the lowest among seven models. The recall rates of the DL, EB, HO, and RML random-effects models are relatively close. Under the third hypothesis, the recall rates of the DSLE2 meta-analysis model are the highest among the seven random-effects models, and the recall rate of the SJ model is the lowest among the seven meta-analysis models. The recall rates of DSLR2, DL, EB, HO, and RML random-effects models are relatively close.

Fig. 7

Plot of recall curves for 7 random-effects meta-analysis models under the first hypothesis

SAR

SAR combines accuracy, AUC (the area under the receiver operating characteristic curve) value, and root mean square error. The SAR value is more robust than a single indicator. We computed the SAR values of seven meta-analysis models under three hypothesis (see Fig. 8, Figure S13 and Figure S14). Under the first hypothesis, the SAR values of the DSLE2 random-effects model are the highest among the SAR values of seven meta-analysis models, and the SAR values of SJ model are the lowest among the seven random-effects models. The SAR values of DL, EB, HO, and RML meta-analysis models are relatively close. Under the second hypothesis, the SAR values of DSLR2 random-effects model are the highest among the SAR values of seven meta-analysis models, and the SAR values of the DSLE2 model are lowest among the seven random-effects models. The SAR values of the DL, EB, HO, RML, and SJ random-effects meta-analysis models are relatively close. The SAR values of the seven meta-analysis models all increase with the increase of the number of studies. Under the third hypothesis, the SAR values of the DSLR2 meta-analysis model are the highest among the SAR values of seven random-effect models. Moreover, the SAR values of the DL, EB, HO, and RML meta-analysis models are relatively close.

Fig. 8

Plot of SAR-value curves for 7 random-effects meta-analysis models under the first hypothesis

Precision-recall F measure

Precision and recall indicators sometimes conflict with each other. In this case, precision and recall indicators need to be considered comprehensively. The most common method is to calculate the F value of precision and recall. We computed the $F$ values of 7 random-effects models under three hypothesis (see Fig. 9, Figure S15 and Figure S16). Under the first hypothesis, the F values of the DSLE2 random-effects model are the highest among the F values of seven random-effects models, and the F values of the SJ meta-analysis model are the lowest among the seven meta-analysis models. The F values of DSLR2, DL, EB, HO, and RML meta-analysis models are relatively close. Under the second hypothesis, the F values of the SJ meta-analysis model are the highest among the F values of seven random-effects models, and the F values of the DSLE2 random-effects model are the lowest among the seven meta-analysis models. The F values of DSLR2, DL, EB, HO, and RML random-effects models are relatively close. Under the third hypothesis, the F values of the DSLE2 random-effects model is the lowest among the F values of seven meta-analysis models.

Fig. 9

Plot of F-value curves for 7 random-effects meta-analysis models under the first hypothesis

Application of DSLE2 random-effects meta-analysis model to lung cancer methylation data

Lung cancer is the second most common cancer worldwide and the most common cancer among men. According to the World Cancer Statistics Center, more than 2.2 million cases of lung cancer occur every year [17]. We collected three sets of lung cancer methylation data in the GEO database: GSE63704, GSE83842, and GSE85845. The EPIC and 450 K methylation data of lung cancer in TCGA database were also collected. A total of 927 samples and 156,680 methylation sites were analyzed. The distribution of the number of significantly differential methylation sites determined by 7 random-effects models is shown in the Fig. 10. The number of significantly differential methylation sites determined by the DSLE2 random-effects model is more than the number of significantly differential methylation sites determined by DSLR2 and is less than the number of significantly differential methylation sites determined by DL, EB, HO, RML, and SJ meta-analysis models.

Fig. 10

Histogram of the number distribution of significantly differential methylation sites identified by 7 random-effects meta-analysis models

The DSLE2 random-effects model identified 59,146 significantly differential methylation sites, of which 2,230 significantly differential methylation sites distributed in 1stExon of 3,039 genes, and 2,003 significantly differential methylation sites were distributed in 3’UTR of 1,660 genes, and 4899 significantly differentially methylated sites distributed in the 5’UTR of 2754 genes, and 21,029 significantly differentially methylated sites distributed in the Body of 7056 genes, and 7500 significantly differentially methylated sites distributed in the TSS1500 of 5113 genes, and 6210 significantly differentially methylated sites distributed in TSS200 of 3963 genes. Research shows that lung cancer is closely related to A2BP1 [8], AACS [18], DNAH10 [19], PINK1 [20], and other genes with significantly differential methylation sites. The significantly differentially methylated sites identified by the DSLE2 meta-analysis model may affect the expression of the corresponding genes.

Application of DSLE2 meta-analysis model to Parkinson’s methylation data

Parkinson’s disease (PD) is a chronic neurodegenerative disease [21]. According to the World Health Organization (WHO), there were 2.5 million and 6.1 million cases of PD worldwide in 1990 and 2016 [22]. However, the number of PD patients increased significantly to 8.5 million, and it is estimated that by 2040, the number of PD patients worldwide may exceed 17 million [23]. In 2019, PD caused 5.8 million disabilities, an 81% increase since 2000, and 329,000 deaths, an increase of more than 100 percent since 2000 [23]. PD causes serious trouble to people’s life, but there is no clear cause and effective treatment. DSLE2 random-effects model was used to further study the causes of PD.

We collected five sets of PD methylation data in the GEO database: GSE72774, GSE72776, GSE111629, GSE145361 and GSE165081. A total of 3,080 samples and 161,261 methylation sites were analyzed. The DSLE2 model identified 26,244 significantly differentially methylation sites (Fig. 11). Most of significantly differentially methylation sites identified by DSLE2, SJ, HO, DSL, DSLR2 meta-analysis methods were same. DSLE2 model independently identified 2181 significantly differentially methylation sites.

Fig. 11

Venn of significantly differential methylation sites identified by 5 random-effects meta-analysis models

We further analyzed the location of significantly differential methylation sites on genes (Figs. 12 and 13). Most of the significant differential methylation sites were located on the gene body, followed by TSS1500 and 5’UTR (Fig. 12). Relatively few Significantly differential methylation sites were located at 3’UTR and the first exon (Fig. 12). Most of the significantly differential methylation sites are uniquely located on a single gene and only a small number of genes contain multiple significantly differential methylation sites (Fig. 13). There are 51, 62 and 122 significantly differentially methylation sites for BCOR, HDAC4 and PRDM16, respectively. Research shows that PD is closely related to BCOR [24], HDAC4 [25], PRDM16 [26], and other genes with significantly differential methylation sites.

Fig. 12

Histogram of significantly differential methylation sites on genes. The horizontal axis refers features, and the vertical axis is the number of significantly differential methylation sites

Fig. 13

Histogram of significantly differential methylation sites on genes. The horizontal axis is the number of significantly differential methylation sites, and the vertical axis is the number of genes

Discussion and conclusion

For high-throughput methylation sequencing data, we proposed the DSLE2 random effects meta-analysis model based on the between-study variance estimator $E_{m}$. Under the alternative hypothesis that the effect sizes of all studies are not 0, the DSLE2 meta-analysis model performs better than the other 6 random-effects models in terms of accuracy, negative prediction rate, recall rate, Matthews correlation coefficient, PCMiss value, and SAR value, F value.

We applied the DSLE2 meta-analysis model to the lung cancer methylation datasets. 59,146 significantly differential methylation sites were identified, with 2,230, 2,230 2003, 4899, 21,029, 7500, and 6210 significantly differentially methylated sites were located on 1stExon, 3'UTR, 5'UTR, Body, TSS1500, and TSS200 of 3039, 1660, 2754, 7056, 5113, and 3963 genes, respectively. Studies have shown that lung cancer is closely related to genes such as A2BP1, AACS, and DNAH10 with significantly differential methylation sites.

We further applied the DSLE2 meta-analysis model to the PD methylation datasets. 26,244 significantly differentially methylation sites, with 952, 776, 2466, 7604, 3897 and 2052 significantly differentially methylated sites were located on 1stExon, 3'UTR, 5'UTR, Body, TSS1500, and TSS200 of 703, 700, 1240, 3059, 2495 and 1267 genes, respectively. Research shows that PD is closely related to BCOR, HDAC4, PRDM16, and other genes with significantly differential methylation sites.

This work presents a new between-study variance estimator for meta-analysis model. First, the primary purpose of the between-study variance estimator is to measure the extent of variability in effect sizes across different studies that is attributable to true differences rather than within-study sampling error. Accurately estimating between-study variance allows us to distinguish between variability due to real differences in study effects and variability due to random error, providing a clearer understanding of the consistency of the effect across studies. Second, meta-analyses synthesize data from multiple studies to derive conclusions with greater statistical power. By accurately accounting for between-study variance, the proposed model enhances the validity of meta-analytic results, leading to more trustworthy conclusions. And it can improve the quality of evidence synthesis, promoting better decision-making in clinical practice, policy, and further research. Third, understanding heterogeneity is key to identifying which treatments work best for specific populations in medical research. A better estimation of between-study variance can help in understanding consistent treatment effects, aiding in the development of personalized therapeutic strategies. Moreover, efficiently allocating research resources requires an understanding of where variability lies. By identifying true sources of heterogeneity, researchers and funding bodies can focus efforts on areas with the most significant impact, optimizing the use of limited resources.

In addition, after significant differential methylation sites identified by DSLE2 meta-analysis model, several software tools and platforms can be used for downstream analysis. These tools help in various aspects such as functional annotation, pathway enrichment analysis and so on. We can use Chip Analysis Methylation Pipeline (ChAMP) to annotate significantly differential methylated regions (DMRs), perform gene ontology (GO) analysis, and integrating with other epigenetic data. GREAT (Genomic Regions Enrichment of Annotations Tool) can be used to annotate and analyze the functional significance of sets of genomic regions, including DMRs, by associating them with nearby genes and functional terms. We can also use GSEA (Gene Set Enrichment Analysis) to do pathway enrichment analysis based on significant differential methylation sites identified by DSLE2.

Data availability

All data used in this paper are available from the GEO database (http://www.ncbi.nlm.gov/geo) and TCGA database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga). The accession numbers of lung cancer methylation data used in this paper are GSE63704, GSE83842, and GSE85845. The EPIC and 450 K methylation data of lung cancer in TCGA database are collected from https://portal.gdc.cancer.gov/analysis_page?app=Downloads. The accession numbers of PD methylation data in the GEO database are GSE72774, GSE72776, GSE111629, GSE145361 and GSE165081.