Quantitative reproducibility analysis for identifying reproducible targets from high-throughput experiments

Background

High-throughput assays are widely used in biological research to select potential targets. One single high-throughput experiment can efficiently study a large number of candidates simultaneously, but is subject to substantial variability. Therefore it is scientifically important to performance quantitative reproducibility analysis to identify reproducible targets with consistent and significant signals across replicate experiments. A few methods exist, but all have limitations.

Methods

In this paper, we propose a new method for identifying reproducible targets. Considering a Bayesian hierarchical model, we show that the test statistics from replicate experiments follow a mixture of multivariate Gaussian distributions, with the one component with zero-mean representing the irreproducible targets.

Results

A target is thus classified as reproducible or irreproducible based on its posterior probability belonging to the reproducible components. We study the performance of our proposed method using simulations and a real data example.

Conclusion

The proposed method is shown to have favorable performance in identifying reproducible targets compared to other methods.

In biological research, high-throughput assays, such as microarrays, are widely used to effectively select potential targets by studying a large number of candidates in a single experiment. However a high-throughput assay is often subject to substantial variability. Reproducibility of high-throughput assays, such as the level of agreement across replicate samples, test sites or data analytical platforms, is a concerned topic in scientific applications, and has been discussed in [1] for microarray and [2] for ChIP-seq technology. Therefore quantitative analysis for the reproducibility of high-throughput assays is an important exercise for evaluating the reliability and robustness of scientific discoveries across studies.

Reproducibility is nonstandard and unsettled across the sciences. Goodman et al. [3] provides a survey on the papers with the word reproducibility included in titles, abstracts and keywords, and concludes that the interpretation of reproducibility varies among different papers. Goodman et al. [3] further allies the word reproducibility in the papers and classifies them into three terms: methods reproducibility, results reproducibility and inferential reproducibility. In [3], methods reproducibility refers to the provision of enough detail about study procedures and data so the same procedures could, in theory or in actuality, be exactly repeated, such as [1] and [2]; results reproducibility refers to obtaining the same results from the conduct of an independent study whose procedures are as closely matched to the original experiment as possible, such as [4] and [5]; Inferential reproducibility refers to the drawing of qualitatively similar conclusions from either an independent replication of a study or a reanalysis of the original study, such as [1] and [2].

In this paper, our reproducibility analysis aims to identify reproducible targets with consistent and significant signals across replicate studies, which belongs to the category of inferential reproducibility as defined in [3]. Our reproducibility analysis is different from meta-analysis, such as [6] and [7]. Meta-analysis combines the data from multiple studies to gain extra power for identifying targets with signals. The identified targets may not necessarily be significant across all studies.

A few methods have been developed for our reproducibility analysis. Hong et al. [8] proposed a permutation based method through estimating the empirical distribution of the rank product. Benjamini & Heller [9] developed a framework for testing partial conjunction hypothesis that the discovery is true in at least u studies out of total n studies. Most recently, [10] proposed a copula mixture model for estimating the irreproducible discovery rate across studies.

However all existing methods potentially have limitations. The permutation based method [8] can be computationally expensive when dealing with a large number of candidates. Benjamini & Heller method [9] aims at identifying candidates with reproduced signals in a few but not all the studies, which is a related but generally weaker goal than ours. The special case of Benjamini & Heller method testing whether signals are reproduced in all studies is identical to using the largest p-value. The copula mixture [10] method builds the copula mixture using the rank transformation of the original data, which might be less powerful than modeling the original data with a proper probabilistic model as in our proposed method. A major drawback of both Benjamini & Heller method [9] and the copula mixture [10] method is that they both use the significant score of signals, such as p-value, without taking into account the directionality of signals, hence is prune to selecting candidates with significant scores but different directions across studies. For example, in the context of two replicate microarray studies with a treatment and a control group, consider genes with significant p-values in both experiments, but are up-regulated in one study and down-regulated in the other. Although those genes have inconsistent signals across studies, both methods will likely classify them as reproducible based on p-values alone. In contrast, our proposed method models the test statistics directly and is expected to correctly classify those genes as irreproducible most of the time.

In this paper, we propose a Bayesian hierarchical model and show the test statistics from replicate studies can be approximated by a mixture of multivariate Gaussian distributions. The proposed Gaussian mixture model classifies the signals into three components: one irreproducible component and two reproducible components for consistent up-regulated and down-regulated signals respectively. The posterior probability of belonging to the reproducible components is used as a measure for reproducibility.

For simplicity, we will introduce our method in the context of microarray studies but it can be generalized to studies of other high-throughput assays. We consider I replicate microarray studies for p genes. In this paper, we focus on the situation of two replicate studies I=2, although our method can be readily extended to the case with more than two studies. We assume a study includes two groups, e.g., the treatment and control group, with sample size equal to n _ik for group k, k=1,2, in the i-th study. Let x _gijk be the normalized and transformed measurement of gene expression of the jth sample from group k for gene g in the i-th study. The test statistics of two-sample unpaired t-test for gene g in the i-th study is

We present an empirical Bayesian hierarchical model to account for various sources of variability. When the sample size n _ik is reasonably large, say n _i1+n _i2≥30, the test statistics d _gi is well approximated by a normal distribution:

where μ _gi is the expected group mean difference for gene g in the i-th study, and \(\delta _{S_{i}}=\tilde {\sigma }_{i}^{-1}(1/n_{i1}+1/n_{i2})^{-1/2}\) with \(\tilde {\sigma }_{i}\) being the common standard deviation for {x _{g i j1}}, j=1,2,…,n _i1 and {x _{g i j2}}, j=1,2,…,n _i2. When the sample size is small, the same procedure as in [11] can be applied to construct z-tests based on two-sample t-tests. For simplicity we assume the within-group between-sample standard deviation is the same for all the genes. The general case can be derived in a similar fashion but a bit more tedious.

For the expected group mean difference μ _gi , we assume it follows

where μ _g is the “true" group mean difference for gene g across all studies and \(\sigma _{g}^{2}\) models the between-study variability due to various experiment conditions.

Furthermore we assume μ _g is from a mixture distribution

where π _i ≥0, i=0,1,2, with π ₀+π ₁+π ₂=1, \(\mu _{G_{1}}>0\) and \(\mu _{G_{2}}<0\). The distribution has three components: the null case where there is no differentially expressed gene, the “up-regulated” case where the treatment stimulates the gene expression, and the “down-regulated” case where the treatment suppresses the gene expression. Generally for microarray studies π ₀≃1. Similar mixture models have been considered in [11–16]. Particularly we choose to model the cluster of up-regulated (or down-regulated) genes with a Gaussian distribution for the computational convenience, same as in [12]. Alternative choices include the semiparametric mixture model in [11, 14], mixture of Gaussian distributions in [13, 15] and mixture of t-distributions in [16].

We can show that the test statistics (d _g1,d _g2) follow a Gaussian mixture model. The derivations are standard by repeatedly applying the law of total expectation and the law of total variance and thus omitted. The mixture model is

where \(\mathcal {N}(\mathcal {\mu }_{l},\Sigma _{l})\) (l=0,1,2) is the biviariate normal distribution with mean vector \(\mathcal {\mu }_{l}\) and covariance matrix Σ _l . Let I ₂ and J ₂ be the identity matrix and the square matrix of ones respectively, both with order 2. This mixture model classify the candidates into three components: \(\mathcal {N}(\mathcal {\mu }_{0},\Sigma _{0})\) is the irreproducible component with zero-mean \(\mathcal {\mu }_{0}=(0,0)^{T}\) and covariance structure \(\Sigma _{0}=\left (\sigma _{g}^{2}+1\right)I_{2}\); \(\mathcal {N}(\mathcal {\mu }_{1},\Sigma _{1})\) and \(\mathcal {N}(\mathcal {\mu }_{2},\Sigma _{2})\) are two reproducible components with \(\mathcal {\mu }_{1}=(\delta _{S_{1}}\mu _{G_{1}}, \delta _{S_{2}}\mu _{G_{1}})>0\) and \(\Sigma _{1}=\left (\sigma _{g}^{2}+1\right)I_{2} + \sigma _{G_{1}}^{2}J_{2}\) representing the up-regulated genes, and \(\mathcal {\mu }_{2}=(\delta _{S_{1}}\mu _{G_{2}}, \delta _{S_{2}}\mu _{G_{2}})<0\) and \(\Sigma _{2}=\left (\sigma _{g}^{2}+1\right)I_{2} + \sigma _{G_{2}}^{2}J_{2}\) representing the down-regulated genes, where the inequalities are meant to be interpreted component-wise.

Note with increased sample sizes or decreased within-group between-sample variability, the mean \(\mathcal {\mu }_{1}\) and \(\mathcal {\mu }_{2}\) of the reproducible components move further away from the origin, making the three components more separable. Also note the test statistics from replicate studies have zero correlations in the irreproducible components; in the reproducible components, the correlations become larger when the between-study variability becomes smaller; for all components, the variance is smaller with less between-study variability, resulting in more separable components.

Under the Gaussian mixture model, the posterior probability of (d _g1,d _g2) belonging to a component is

where ϕ(·|·) is the density function of bivariate normal distribution. According to [10], the posterior probability of being in the irreproducible/null component p _i0 can be introduced as the individual significant score, namely local false discovery rate. When p _g0 is less than a significant level α, gene g is classified as reproducible.

Next, we consider estimation of the unknown parameters

in the mixture model (4) to get the estimate of p _g0 for individual genes. It is natural to use the expectation-maximization (EM) algorithm to estimate \(\mathcal {\theta }\) by maximizing the log-likelihood of the data [17], i.e.,

In our algorithm, we start with some initials value for the parameters \(\mathcal {\theta }^{0}\), then iterate between two steps: (1) Evaluate the current posterior probabilities p _gl using the current parameters; (2) Maximize the likelihood estimator given current posterior probabilities. The details of the EM procedures are provided in Appendix. Multiple random initial vaues are used to avoid being trapped at the local maximum.

In this section, we present numerical simulations to illustrate the performance of our proposed method compared to three existing methods, the copula mixture model [10], Benjamini & Heller method [9], and the rank product method [8]. We use the following model to simulate data

From this model, the mean expression level of gene g for group 1 of study s is modeled as μ _{g s1}=μ+α _g +β _i +(α β) _gi , where μ is the overall mean; α _g is the main effect of gene g; β _i is the main effect of study i; (α β) _gi is the gene-study interaction. We set μ=0, \(\alpha _{g} \sim \mathcal {N}(0,1)\), β _i =0.1, and \((\alpha \beta)_{gi} \sim \mathcal {N}(0,0.5^{2})\). For non-differentially expressed genes, the mean expression level for both groups are the same, i.e., μ _{g s1}=μ _{g s2}. For differentially expressed genes, (8) models the difference between the two comparison groups as μ _{g i2}−μ _{g i1}=δ+γ _g +(γ β) _gi , where δ is the fixed effect of group difference; γ _g is the effect of gene on the group difference; (γ β) _gi is the gene-study interaction of the group difference. We set δ=0, generate γ _g from \(\mathcal {N}(2, 0.5^{2})\) or \(\mathcal {N}(-2, 0.5^{2})\) to mimic two possible directions of signals, \((\gamma \beta)_{gi}\sim \mathcal {N}(0,0.5^{2})\). ε _gijk is the random error term, and following the distribution \(\mathcal {N}(0,0.5^{2})\).

For each simulation run, we generate 2 studies. Each study has two groups with 10 samples per group. We generate G=5000 genes per sample and choose the proportions of reproducible genes (γ) from (80%, 60%, 40%, 20%, 10%, 5%, 1%). We apply the proposed method and the three existing methods to the simulated data, and classify the genes as reproducible based on two commonly used significant levels (α) 0.05 and 0.1. The performance of the four compared methods is evaluated by three criteria, i.e., sensitivity, specificity and misclassification rate. Results from 50 simulations are summarized in Tables 1, 2 and 3 respectively. The results shows our proposed method performs the best among the four methods with the smallest misclassification rates (Table 1), highest sensitivity (Table 2) and highest specificities (Table 3).

In this section, we illustrate our proposed method using a real example. This example includes two microarray studies [18] and [19] comparing idiopathic pulmonary fibrosis (IPF) samples with healthy control samples. Data from both studies are obtained from Gene Expression Omnibus [20]. GSE 28042 [18] measures profiles of peripheral blood mononuclear cell (PBMC) for 75 IPF samples and 16 control samples through GeneChip Human 1.0 exon ST arrays, and GSE 33566 [19] measures profiles of peripheral blood RNA for 93 IPF patients and 30 control samples through Agilent Whole Human Genome Oligonucleotide Microarrays. We only consider the overlap 17708 common genes for reproducibility analysis.

We apply our proposed method, the copula mixture model [10] and Benjamini & Heller method [9]. The rank product method [8] is too computationally intensive to be applied to this example and thus excluded from this study. Figures 1, 2 and 3 show the results of selected reproducible genes from the three compared methods respectively (green). In all three figures, the x axis represents the test statistics from GSE 28042 [18], and the y axis represents the test statistics from GSE 33566 [19]. The top 500 reproducible genes selected by three methods are highlighted in green. As shown in Fig. 1, our proposed method only selects genes with consistently significant signals in both studies. Benjamini & Heller method [9] incorrectly identifies 23 genes (the upper left and bottom right corners of Fig. 2) as reproducible, which actually have opposite directions in two studies. The complete list of the 23 genes incorrectly selected by Benjamini & Heller method [9] is provided in Table 4. The copula mixture model [10] selects 7 genes (Table 5) with opposite directions of signals. It’s also noted that the copula mixture model [10] appears to be less powerful in separating the irreproducible and reproducible genes and has incorrectly selected some insignificant genes (see the center of Fig. 3), likely resulting from the rank transformation. Overall, our method performs favorably in identifying reproducible genes.

Bivariate plot of test statistics from two studies. The x axis represents the test statistics from GSE 28042 study [18], and the y axis represents the test statistics from GSE 33566 [19]. The green points are the top 500 reproducible genes selected by the proposed method

Full size image

Bivariate plot of test statistics from two studies. The x axis represents the test statistics from GSE 28042 study [18], and the y axis represents the test statistics from GSE 33566 [19]. The green points are the top 500 reproducible genes selected by the copula mixture model [10]

Full size image

Bivariate plot of test statistics from two studies. The x axis represents the t-statistics from GSE 28042 study [18], and the y axis represents t-statistics from GSE 33566 [19]. The green points are the top 500 reproducible genes selected by Benjamini & Heller method [9]

Full size image

This paper proposes a new method for identifying consistent and significant signals across replicate high-throughput experiments. Existing methods ignore the directionality of signals, and can incorrectly identify signals with opposite directions as reproducible ones. Our proposed method considers both the significant scores and directions of signals by modeling the test statistics directly, leading to improved performance in selecting reproducible candidates. When the proposed method is applied to a real data example for identifying reproducible genes in studies of idiopathic pulmonary fibrosis samples, it is shown to have better performance in detecting significant and reproducible genes compared to other methods. Simulations also demonstrate that our method compares favorably to the existing methods.

Expectation-maximization (EM) algorithm to estimate model parameters

The algorithm for estimating \(\mathcal {\theta }\) in (6) is an iterative algorithm between Expectation steps and maximization step. We use \(\widehat {\mathcal {\theta }}^{v}\) to denote the estimate at vth iteration. The algorithm includes the following steps:

• Step 1: Initial Values Generate the initial values for \(\mathcal {\theta }\) and denote it as \(\widehat {\mathcal {\theta }}^{0}\)

• Step 2: Expectation-Step Continue from the vth iteration step with the estimate \(\widehat {\mathcal {\theta }}^{v}\). We can obtain the estimated posterior probability \(\widehat {p_{gl}}^{v}\) of (d _g1,d _g2) from (5) by

  $$\begin{array}{@{}rcl@{}} \widehat{p_{gl}}^{v}&=&\frac{\widehat{\pi_{l}}^{v}\phi\left(d_{g1},d_{g2}|\widehat{\mathcal{\mu}_{l}}^{v},\widehat{\Sigma_{l}}^{v}\right)}{\sum_{\ell=0,1,2}\widehat{\pi_{\ell}}^{v}\phi\left(d_{g1},d_{g2}|\widehat{\mathcal{\mu}_{\ell}}^{v},\widehat{\Sigma_{\ell}}^{v}\right)},l=0,1,2. \end{array} $$

  (9)

• Step 3: Maximization-Step Update the parameter \(\widehat {\mathcal {\theta }}^{v+1}\) by maximizing the log-likelihood function \(\ell (\mathcal {\theta })\) in (7) given the current estimated posterior probability \(\widehat {p_{gl}}^{v}\). The estimated parameters from the maximization are

  $$\begin{aligned} \widehat{\pi_{l}}^{v+1}&=\sum_{g=1}^{p}\widehat{p_{gl}}^{v}/p,l=0,1,2.\\ \widehat{\mu_{1}}^{v+1}&=\left(\widehat{\mu_{11}}^{v+1},\widehat{\mu_{12}}^{v+1}\right)= \left(\frac{\sum_{g=1}^{p} \widehat{p_{g2}}^{v} d_{g1}}{\sum_{g=1}^{p} \widehat{p_{g2}}^{v}}, \frac{\sum_{g=1}^{p} \widehat{p_{g2}}^{v} d_{g2}}{\sum_{g=1}^{p} \widehat{p_{g2}}^{v}}\right)\\ \widehat{\mu_{2}}^{v+1}&=\left(\widehat{\mu_{21}}^{v+1},\widehat{\mu_{22}}^{v+1}\right)= \left(\frac{\sum_{g=1}^{p} \widehat{p_{g3}}^{v} d_{g1}}{\sum_{g=1}^{p} \widehat{p_{g3}}^{v}}, \frac{\sum_{g=1}^{p} \widehat{p_{g3}}^{v} d_{g2}}{\sum_{g=1}^{p} \widehat{p_{g3}}^{v}}\right) \\ \widehat{\sigma_{g}^{2}}^{v+1}&= \frac{\sum_{g=1}^{p} \widehat{p_{g1}}^{v}\left(d_{g1}^{2}+d_{g2}^{2}\right)}{2\sum_{g=1}^{p} \widehat{p_{g1}}^{v}} -1 \\ \widehat{\sigma_{G_{1}}^{2}}^{v+1}&= \frac{\sum_{g=1}^{p} \left[\widehat{p_{g2}}^{v} \left(d_{g1}-\widehat{\mu_{11}}^{v+1}\right)^{2}+\widehat{p_{g2}}^{v} \left(d_{g2}-\widehat{\mu_{12}}^{v+1}\right)^{2}\right]}{2\sum_{g=1}^{p} \widehat{p_{g2}}^{v}} \\ &\quad-\frac{\sum_{g=1}^{p} \widehat{p_{g1}}^{v}\left(d_{g1}^{2}+d_{g2}^{2}\right)}{2\sum_{g=1}^{p} \widehat{p_{g1}}^{v}} \\ \widehat{\sigma_{G_{2}}^{2}}^{v+1}&= \frac{\sum_{g=1}^{p} \left[\widehat{p_{g3}}^{v} \left(d_{g1}-\widehat{\mu_{21}}^{v+1}\right)^{2}+\widehat{p_{g3}}^{v} \left(d_{g2}-\widehat{\mu_{22}}^{v+1}\right)^{2}\right]}{2\sum_{g=1}^{p} \widehat{p_{g3}}^{v}} \\ &\quad-\frac{\sum_{g=1}^{p} \widehat{p_{g1}}^{v}\left(d_{g1}^{2}+d_{g2}^{2}\right)}{2\sum_{g=1}^{p} \widehat{p_{g1}}^{v}} \end{aligned} $$

• Step 4: Solution The algorithm continues between Expectation-Step and Maximization-Step until the following two conditions are satisfied.

  1. 1.

     The difference between \(\widehat {\mathcal {\theta }}^{v} \) and \(\widehat {\mathcal {\theta }}^{v+1}\) is less than a small value δ ₁ for all their elements;

  2. 2.

     The change in log-likelihood function \(\ell (\mathcal {\theta })\) between two consecutive iterations does not exceed a small value δ ₂.

We would like thank referees for their time on reviewing this manuscript.

Availability of data and materials

All data are from Gene Expression Omnibus [20].

Affiliations

1. Department of Biostatistics, Columbia University, New York, NY, USA

   • Ying Liu

2. Sanofi, Framingham, MA, USA

   • Wenfei Zhang
   • , Mindy Zhang
   • , Cheng Zhu
   •  & Yuefeng Lu

Contributions

All authors equally distributed. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Wenfei Zhang.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Reprints and Permissions