Increasing reproducibility, robustness, and generalizability of biomarker selection from meta-analysis using Bayesian methodology.

A major limitation of gene expression biomarker studies is that they are not reproducible as they simply do not generalize to larger, real-world, heterogeneous populations. Frequentist multi-cohort gene expression meta-analysis has been frequently used as a solution to this problem to identify biomarkers that are truly differentially expressed. However, the frequentist meta-analysis framework has its limitations-it needs at least 4-5 datasets with hundreds of samples, is prone to confounding from outliers and relies on multiple-hypothesis corrected p-values. To address these shortcomings, we have created a Bayesian meta-analysis framework for the analysis of gene expression data. Using real-world data from three different diseases, we show that the Bayesian method is more robust to outliers, creates more informative estimates of between-study heterogeneity, reduces the number of false positive and false negative biomarkers and selects more generalizable biomarkers with less data. We have compared the Bayesian framework to a previously published frequentist framework and have developed a publicly available R package for use.

Introduction

With the advent of high-throughput transcriptomics, researchers have been able to profile gene expression in millions of samples at low costs, opening many new avenues of research. However, single-cohort studies with low-to-moderate statistical power are one of the four horsemen of irreproducibility in biomedical research [1]. Recent studies have estimated that only 10–25% of biomedical studies are reproducible [2-4]. This is especially a problem in biomarker studies, where differentially expressed genes (DEGs) between subjects with disease of interest and controls rarely generalize to the real-world patient populations. One of the reasons for this lack of generalizability is not accounting for heterogeneity representative of the real-world patient population.

Broadly, there are three sources of heterogeneity in the real-world patient population: biological (age, sex, tissue, cell type), clinical (treatment, disease duration, comorbidities), and technical (experimental protocol, batch effects). Traditionally, in a single cohort analysis, these sources of heterogeneity reduce the statistical power, requiring a large number of samples. Therefore, single cohort studies strive to increase statistical power by limiting heterogeneity as much as possible. However, this reduction in heterogeneity in single cohort studies leads to reduced generalizability to heterogeneous, real-world patient populations. We have repeatedly shown that leveraging heterogeneity across independent cohorts using a frequentist meta-analysis approach can identify robust disease signatures that are diagnostic and prognostic, and have been translated into a point-of-care test for clinical use [5,6]. We have previously identified best practices for gene expression meta-analyses, such as the number of studies and samples needed [7].

Despite its repeated success in identifying robust disease signatures, the frequentist approach has its limitations. First, previous work has shown that approximately 4–5 datasets with about 250 samples are needed to perform a successful frequentist meta-analysis [7]. However, several diseases simply do not have enough samples or datasets publicly available for successful integration using this guideline. Second, the statistic used to estimate and summarize effect sizes (e.g., Cohen’s d, Hedge’s g) can be susceptible to outlier samples within a subset of studies, resulting in misleading effect size estimates. Finally, frequentist approaches rely on multiple hypotheses corrected p-values, which are shown to be substantially underestimated [7]. Bayesian meta-analysis approaches have the potential to overcome these limitations. For example, Bayesian estimation has previously been shown to be more outlier resistant than traditional hypothesis testing [8]. Importantly, unlike frequentist meta-analysis, adjusting for multiple comparisons is not required for Bayesian approaches and yields more efficient and reliable estimates of effect [9].

We compared a frequentist approach with a new framework utilizing Bayesian approximation supersedes the t-test (BEST) for the meta-analysis of transcriptome using multiple independent datasets from humans with different diseases [8]. Here, we show that using this Bayesian approach, we are able to: 1) select more generalizable and robust biomarkers with fewer datasets, 2) be robust to outliers, 3) create better estimates of between study heterogeneity for biomarker selection, and 4) reduce the number of false positives and false negative genes for classification. This framework has also been developed into an R package, bayesMetaIntegrator, that is publicly available for use (https://github.com/Khatri-Lab/bayesMetaIntegrator).

Results

Bayesian meta-analysis is resistant to outliers and provides a better estimate of heterogeneity in gene expression meta-analysis

We investigated whether the two meta-analysis approaches, Bayesian ( and frequentist, identified the same or different set of genes using four publicly available asthma bronchial epithelial cell gene expression datasets—differentiating samples from asthma patients and healthy controls () [10-12]. Although the summary effect sizes were highly correlated between approaches (r = 0.94, p<2.2e-16; ), the Bayesian approach consistently estimated higher between-dataset heterogeneity, τ2, than the frequentist approach (). While the frequentist approach found a large number of genes (26%) with no between-dataset heterogeneity, the Bayesian approach did not assign near-zero heterogeneity (τ2 < 0.01) to any genes. This difference in τ2 between the two approaches is due to how it is estimated. In a frequentist meta-analysis, high within-study heterogeneity leads to wider confidence intervals, which in turn drowns out the between-study heterogeneity. In contrast, Bayesian meta-analysis uses a probabilistic distribution to represent τ2 instead of a confidence interval, resulting in more conservative estimates of heterogeneity. The range of between-datasets heterogeneity was substantially higher for the Bayesian approach compared to the frequentist approach, further suggesting a Bayesian approach is more conservative than a frequentist approach.

Bayesian gene expression schematic.

We first use the BEST framework to estimate the posterior distribution of effect size between cases and controls for each gene in each dataset. Then we combine the distributions from independent studies using a gaussian hierarchical model, estimating both the pooled effect size and between-study heterogeneity in the process. Following, we estimate the probability of a gene being upregulated or downregulated based on the pooled effect size distribution.

Asthma Meta-Analysis Comparison.

(A) Comparison of pooled effect sizes by the frequentist and Bayesian meta-analysis for all genes. Pearson correlation (r) = 0.94. (B) Comparison of between-dataset heterogeneity, τ2, by the frequentist and Bayesian meta-analysis for all genes. Pearson correlation (r) = 0.80. (C) Comparison of within-dataset effect size estimates by the frequentist and Bayesian meta-analysis for all genes. (D) Comparison between false discovery rate from the frequentist meta-analysis and probabilities from the Bayesian meta-analysis. (E) The x axes represent standardized mean difference between patients with asthma and healthy controls, computed as Hedges’ g, in log2 scale. The size of the blue rectangles is proportional to the standard error of mean difference in the study. Whiskers represent the 95% confidence interval. The orange diamonds represent overall, combined mean difference (summary effect size) for a given gene. Width of the orange diamonds represents the 95% confidence interval of overall mean difference. (F) Expression of SLC6A4 in GSE64913 and FANCI in GSE18965 in patients with asthma and healthy controls.

Interestingly, despite the high correlation between summary effect sizes, within-dataset effect size correlations ranged widely from 0.87 to 0.95. Closer examination found that both approaches differed substantially in effect size estimates for a subset of genes. For example, in GSE64913, a subset of genes has an effect size of 0 when using the Bayesian approach, but a non-zero effect size when using a frequentist meta-analysis (). Therefore, we investigated whether these differences in effect size estimates and between-dataset heterogeneity will lead to identification of different set of differentially expressed genes (DEGs).

We used false discovery rate (FDR) for the frequentist approach and Bayesian probability for the Bayesian approach as the measures of statistical significance. Low correlation between both measures (r = 0.37, p<2.2e-16) suggested they identified different sets of DEGs. As the number of significant genes by either approach increased, the number of genes identified as significant by both methods also increased (). Out of 1,117 DEGs that were statistically significant by either approach, 545 genes (48.8%) were statistically significant by both approaches (FDR < 10% and Bayesian probability < 0.1). The remaining 572 genes were significant when using either the frequentist approach (235 genes, 21%) or the Bayesian approach (337 genes, 30.2%) ().

Importantly, for genes that were significant by either approach but not both, we found that the Bayesian approach is more robust to outlier samples within a single dataset than the frequentist approach. For example, SLC6A4 was significant by the frequentist approach (FDR<4%), but not by the Bayesian approach (p = 0.24). Comparison of effect sizes from both approaches showed that for GSE64913, the frequentist approach estimated statistically significant non-zero effect size, whereas the Bayesian approach estimated non-significant effect size (). Further analysis showed that although 59 out of 70 samples in GSE64913 had identical expression values for SLC6A4, irrespective of estimate of statistically significant effect size by the frequentist approach was driven by only 2 out of 70 samples (). In contrast, the Bayesian approach correctly estimated near-zero effect size due to its reliance on parameter estimation and sampling, and was not confounded by a small number of outliers. This observation demonstrated that the Bayesian approach reduced false positives.

In contrast, another gene, FANCI, was statistically significant by the Bayesian approach (p = 0.01), but not by the frequentist approach (FDR = 19%). Although both approaches estimated the effect size for FANCI in the same direction, the difference in statistical significance was due to the difference in estimated effect size in a single dataset (GSE18965) (). The frequentist approach estimated negative effect size for FANCI in GSE18965, which was driven by a single healthy control sample (), which in turn led to higher between-study heterogeneity for the gene and the summary effect size for the gene being statistically insignificant. In contrast, the Bayesian approach was not confounded by a single sample, correctly estimated its effect size as positive in GSE18965, and identified it as being statistically significant overall. Collectively, these results show that the Bayesian meta-analysis approach is robust to outliers, which in turn reduces false positives and false negatives.

Comparing Bayesian and frequentist meta-analysis methods

A key difference between the frequentist and Bayesian meta-analysis approaches is how between-study heterogeneity affects estimates of summary effect size and its statistical significance. To investigate the effect of between-study heterogeneity, we simulated data across 5 independent studies (Methods) such that either the effect size, its variance, or both changed for one study (). Increasing the effect size for one study without changing its variance estimate, which simulated adding a study with high certainty of a strong positive effect, statistical significance reduced for the frequentist approach (i.e., higher FDR), but increased for the Bayesian approach (i.e., lower probability), although the between-study heterogeneity increased for both approaches (). When we increased only variance or both effect size and variance, the summary heterogeneity and statistical significance decreased for both approaches (). These results, combined with robustness of the Bayesian approach to outliers within a dataset, are desired characteristics for identifying generalizable signal across heterogeneous datasets.

Simulated Meta-Analysis.

A) Using simulated data in five studies, we explored the effect of changing the effect size and variance of a single study on each model’s pooled effect size, heterogeneity (τ2) and significance. Here, we show the baseline meta-analysis. We picked a scenario where both Bayesian and frequentist are borderline significant. B) Increasing effect size. We increased the effect size of a single study–creating high certainty of a strong positive effect of that study. The Bayesian model became more significant, whereas the frequentist model became less significant. Both increased in effect size and τ2. C) Increasing variability. We increased the variability of a single study. We found that both methods behave similarly–with effect size, significance and τ2 decreasing. D) Increasing variability and effect size. We increased both the variability and effect size of a single study. We found that both methods behave similarly–with effect size, significance and τ2 decreasing.

Comparing the predictive performance of Bayesian and frequentist meta-analysis

Given that the Bayesian approach is robust to outlier samples within a study and better estimates between-study heterogeneity, we investigated whether it would identify DEGs that are more generalizable to unseen data than the frequentist approach. To investigate this, we applied both meta-analysis approaches to transcriptome profiles from patients with cardiomyopathy (14 datasets, 1039 samples) [13-23] or tuberculosis (27 datasets, 3069 samples) [24-45]. For both diseases, we identified most differentially expressed genes by successively increasing the number of datasets and compared their discriminatory power in unused datasets ( and Methods). For both diseases, we randomly selected N datasets 100 times and applied both meta-analysis approaches to each set of randomly selected N datasets. For each iteration, we selected DEGs with absolute summary effect size > 0.6 and used the top 10 genes with the lowest FDR or Bayesian probability. We used difference of geometric means of over- and under-expressed genes as a classifier in unseen datasets to distinguish cases from healthy controls. We chose difference of geometric mean as a classification model because such a classifier has been repeatedly demonstrated to be generalizable and has been translated in a point-of-care test [5,6,46,47]. When we varied the effect size threshold (0.4 to 1.1) and the number of selected genes (10 to 200), while keeping the number of datasets used for analysis constant at 4, the genes selected using the Bayesian approach consistently led to higher AUC ().

Comparing gene signatures from frequentist and Bayesian meta-analysis.

A) Selection method—for both diseases, we randomly selected N datasets 100 times. For each iteration, we selected DEGs with absolute pooled effect size > 0.6 and selected 10 genes with the smallest FDR or Bayesian probability. The signature is then created using a geometric mean signature score and validated on completely left-out datasets. B) Average AUC–the Bayesian method picked gene signatures that are more generalizable—achieving higher AUCs for both diseases on left-out datasets. The Bayesian method also requires less datasets to achieve higher classifier performance. C) Average τ–we find that the between study heterogeneity (τ) is significantly lower for the gene signatures selected by the Bayesian method compared to the frequentist, suggesting the Bayesian method is selecting genes that have lower heterogeneity. Due to the frequentist method calculating τ as 0 for a large percentage of genes, Bayesian tau is displayed for all genes. D) Mean pooled effect size–we do not see consistent trends in the average effect sizes of the gene signatures, except for the tuberculosis analysis where the frequentist method tended to pick higher effect size genes with larger Ns. **** < 0.0001, *** < 0.001, ** < 0.01, * < 0.05, ns: > 0.05.

First, we compared area under the receiver operating characteristic (AUROC) curves in the unused datasets in each iteration as a proxy for identifying generalizable gene signatures for both meta-analysis approaches. For both diseases, irrespective of the number of datasets used, the DEGs identified using the Bayesian approach had consistently higher AUROC in unseen datasets than using the frequentist approach. When using 2 out of 27 datasets to identify DEGs for tuberculosis, there was no significant difference in AUROC between the two approaches, which suggests that N = 2 does not represent the heterogeneity across the other 27 datasets. Interestingly, for both diseases, the median AUROCs for the Bayesian approach using 3 datasets was always equal or greater than the median AUROCs for the frequentist approach using substantially larger number of datasets (9 datasets for cardiomyopathy, 14 datasets for tuberculosis). This result demonstrated that the DEGs identified using the Bayesian approach are more generalizable to previously unseen data than those identified using the frequentist approach. Further, it also suggests that using as few as 3 datasets may be sufficient to identify robust gene signatures of a disease using the Bayesian meta-analysis approach.

Next, we investigated between-study heterogeneity, τ2, and effect sizes for the DEGs identified by the two approaches and whether those differences explained more generalizability for the Bayesian approach. For both diseases, median τ2 decreased with the increasing number of datasets (). The DEGs identified by the Bayesian approach had significantly lower τ2 than those identified by the frequentist approach. Interestingly, for the frequentist approach as the number of datasets increased, variability in τ2 increased for tuberculosis analysis but not cardiomyopathy analysis. However, there was no consistent difference in the summary effect sizes of the genes selected by either approach. For example, for cardiomyopathy datasets, the difference in summary effect sizes of DEGs were mostly non-significant, whereas for tuberculosis datasets, there was no statistically significant difference in summary effect sizes, when using a smaller number of datasets (N≤6), but the frequentist approach selected genes with significantly higher summary effect sizes for a larger number of datasets (N>6). Interestingly, for tuberculosis analysis, higher summary effect size and variability in τ2 with the increased number of datasets suggest that the DEGs may be affected by outlier datasets, which in turn reduces their generalizability to unseen datasets with reduced AUROCs.

To further investigate the differences between gene sets created using Bayesian or frequentist meta-analysis, we performed an overall meta-analysis with all of the Asthma, Cardiomyopathy and Tuberculosis datasets. Calculating the Jaccard similarity of the gene sets, we found that for the smallest meta-analysis, asthma, the Jaccard similarity between the models was low for the top 500 genes, only reaching 0.3. For the larger meta-analyses, Cardiomyopathy and Tuberculosis, we found that the Jaccard similarity increased up to 0.5 (). In each case, the Jaccard similarity increased until reaching a plateau. We also compared the gene sets using pathway analysis using ReactomePA [48] (). The pathways represented were similar in both analysis types for all three analyses, which suggested that although Bayesian meta-analysis identified more generalizable genes that have higher discriminatory power, it is still identifying the same biological pathways as frequentist meta-analysis.

Collectively, our results show that compared to the frequentist approach, the Bayesian approach for meta-analysis identifies genes with lower between-study heterogeneity and comparable summary effect sizes, and is robust to outlier samples, which in turn leads to more generalizable classifier for unseen datasets. Our results also suggest that the Bayesian approach requires lower number of datasets to identify generalizable DEGs compared to the frequentist approach.

Discussion

We performed three gene expression meta-analyses to compare the Bayesian and frequentist meta-analysis approaches. Using dozens of publicly available gene expression studies, we found that Bayesian approach tends to identify differentially expressed genes that have lower between-dataset heterogeneity and higher discriminatory power, which leads to more generalizable classifiers. Importantly, we found that the Bayesian approach consistently required lower number of datasets than the frequentist approach.

Several factors contribute to drive these effects for the Bayesian meta-analysis approach. First, our analysis showed that the Bayesian approach is resistant to outliers due to the t-distribution underlying the estimation of the effect size distribution per dataset. Second, the Bayesian approach uses probabilistic distributions to represent effect size in each dataset as opposed to confidence intervals that the frequentist approach relies on. When the number of datasets used for meta-analysis is small and the within-dataset variability is higher, confidence intervals tend to be wider, leading the frequentist approach to estimate the between-study heterogeneity as zero or very low for a large number of genes the large within dataset heterogeneity. In contrast, the use of probabilistic distribution leads to conservative estimates of between-study heterogeneity. For biomarker discovery, this is preferable, for we seek to find the biomarkers that have the smallest between-study heterogeneity across all datasets in our analysis. Finally, the p-values in the frequentist approach represent the probability of observing data under a hypothesis of no effect and must be multiple hypothesis adjusted, whereas Bayesian probability represents the posterior belief of the difference between groups and require no multiple hypothesis adjustment [49].

Although the Bayesian meta-analysis tends to perform better than the frequentist method at finding consistently differentially expressed genes across studies with low heterogeneity, one area in which it would be less advantageous is the unsupervised identification of subgroups within patient populations. In this case, we would want to select for genes that have high effect sizes when compared to controls, but also have heterogeneity to separate between cases. For this task, the frequentist method would likely be more effective than the Bayesian method described in this study. However, priors and probability calculations can be adjusted, providing Bayesian meta-analysis the flexibility to succeed in many different scenarios. A limitation of our study is that we used minimally informative priors for all estimation in order to produce the most accurate estimates of effect size and heterogeneity. However, depending on the context, these priors could be changed. For example, in the case of finding biomarkers for diagnostic use, the prior for τ2 could be changed from a uniform to a monotonically increasing function. This would in turn create a form of regularization, pushing τ2 to be estimated as larger and increasing the effect of heterogeneity on Bayesian probability estimate. For subgroup identification and clustering of patients, one could similarly adjust priors and probability estimates to select genes that have high effect size and moderate-to-high within and/or between study heterogeneity. This shows the true flexibility and potential adaptability of the Bayesian framework for different uses.

Methods

Dataset selection

We used publicly available transcriptome data from the NCBI GEO for three diseases: (1) 223 samples across 4 datasets from healthy controls and patients with asthma () [10-12], (2) 1039 samples across 14 datasets from healthy controls and patients with cardiomyopathy () [13-23], and (3) 3069 samples across 27 datasets from healthy controls and patients with tuberculosis (TB) () [24-45]. Each dataset was appropriately normalized and log2 transformed, if not already in log scale. We removed genes that were not present in at least half of the datasets for a given disease.

Frequentist meta-analysis

We used the frequentist meta-analysis implemented in MetaIntegrator, which uses random effects inverse variant model, for comparison with the Bayesian method [50]. Briefly, MetaIntegrator computes a Hedge’s g as an effect size for each gene in each dataset. The effect sizes are combined using random effects inverse variance model the DerSimonian-Laird method, and the corresponding p-value is estimated using a standard normal distribution, which is corrected for multiple hypotheses testing using the Benjamini-Hochberg FDR adjustment [51]. Following these calculations, the top genes are selected by using FDR and effect size thresholds. We used difference between geometric mean of over-expressed genes to that of under-expressed genes as a classifier because it has been repeatedly shown to be more generalizable across datasets [52] and has also been translated in a point-of-care test [5,6,46,47].

Bayesian meta-analysis–dataset effect size calculation

The first step of the Bayesian meta-analysis pipeline involves creating an effect size distribution for each case and control for each gene in each dataset (). We used the BEST [8] framework with default parameters and priors for this purpose. The BEST framework estimates the posterior distribution of effect size between cases and controls for each gene in each dataset by assuming that the data is independently distributed and comes from a t distribution with different mean (μ) and standard deviation (σ) parameters for each group. Then, we combine the distributions from independent studies using a gaussian hierarchical model, estimating both the pooled effect size and between-study heterogeneity in the process. Overall normality parameter (v) that denotes the size of tails of the t distribution and the level of normality. Overall, the BEST framework estimates 5 parameters: μ1, μ2, σ1, σ2 and v using minimally informative priors. μ1 and μ2 are the population means of cases and controls and are parameterized with a wide normal prior with a large standard deviation. σ1 and σ2 are the population standard deviations of cases and controls and are parameterized with a broad uniform. The normality parameter, v, has a broad, shifted exponential prior. Following the parameter estimation, we calculated the effect size as a standardized mean difference:

For this study, we ran all individual datasets with 2000 steps with 400 for model burn-in with 3 chains to ensure convergence, which is then calculated using [53]. We removed any genes with an greater than 1.1 from the dataset.

Bayesian meta-analysis–pooling step

Following the dataset effect size distribution estimation, a pooling step is performed to estimate an overall pooled effect size using a hierarchical model. The effect size distribution from each dataset for each gene is assumed to be normally distributed with mean μi and σi. To calculate the pooled effect size, we use each of the calculated dataset effect size distributions and assume that they are sampled from an overall, pooled normal distribution represented as . Both μ and τ are parameterized with minimally informative priors: μpooled as and τ as and parameters are estimated using Gibbs sampling [54]. We chose priors for effect size and between-study heterogeneity using 122 previous gene expression meta-analyses ( and ) [50]. Sensitivity analysis of the priors found that when we varied the parameters of the Normal priors for effect size, the Bayesian probabilities for differential expression remained concordant ( and ). For between-study heterogeneity, we found that using Uniform(0,1) as prior had lower Bayesian probability estimates (), whereas using Uniform(0,2) or Uniform(0,3) as prior had highly concordant Bayesian probability estimates (). Hence, the final gene rankings and posterior probabilities did not change by widening the priors further.

For this study, all pooling steps were run with 5000 steps with 1000 for burn-in with 3 chains. The convergence parameter () is calculated using the chains and any genes with an greater than 1.1 are removed. Using the hierarchical model structure, once the pooled distribution is estimated, we adjusted the individual dataset effect size distributions based on the summary distribution, akin to a random effects model in frequentist meta-analysis. To calculate statistical significance, we calculate the probability that a certain gene is upregulated or downregulated, calculating Pr(μpooled < 0) or Pr(μpooled > 0), respectively. This is done with a standard cumulative density function for a normal distribution.

Simulated data for comparison of frequentist and Bayesian meta-analysis

Random study data was used to compare the Bayesian and frequentist methods. The frequentist random-effects meta-analysis was run using the metagen function from the meta package [55] in the R using default parameters. The Bayesian model was run as with the parameters described above, with 50000 steps with 10000 for burn-in.

Comparison of AUCs, effect sizes, and heterogeneity between frequentist and Bayesian meta-analysis

To compare the methods, we used the Tuberculosis and Cardiomyopathy cohorts, as defined above. We removed genes that were not present in at least half of the datasets. For both diseases, we randomly selected N datasets 100 times (or 91 times in the case of Cardiomyopathy: N = 2) and applied both meta-analysis approaches Bayesian (bayesMetaIntegrator) and frequentist (MetaIntegrator) to each set of randomly selected N datasets. For each iteration, we filtered to differentially expressed genes with absolute summary effect size > 0.6 and used the top 10 genes with the smallest FDR or Bayesian probability. We used difference of geometric means of over- and under-expressed genes as a classifier in unseen datasets to distinguish cases from healthy controls, as defined below. For the tau comparison, we report the Bayesian estimate of tau for each gene due to frequentist meta-analysis reporting large numbers of genes as 0. For pathway analysis, the ReactomePA package [48] was used to perform pathway analysis on all three diseases, using the most significant genes for both methods– 1000 genes for Cardiomyopathy and Tuberculosis and 500 genes for Asthma. A p-value cutoff of 0.2 was used for Cardiomyopathy, 0.2 for Asthma and 0.05 for Tuberculosis.

A) Using 122 previous gene expression meta-analyses we observed the pooled effect size and tau to determine our initial priors. Our priors, picked to be minimally informative, are shown in red. B) Sensitivity analysis of effect size prior using Asthma data. C) Sensitivity analysis of tau prior using Asthma data.

(EPS)

Click here for additional data file.

A) Using N = 4 datasets and top 10 statistically significant genes, we examined the effect of effect size thresholds using the AUC performance on all other left-out datasets. We find that the Bayesian model consistently outperforms the frequentist at effect sizes < 1. B) Using N = 4 datasets and effect size threshold of .6, we examined how the number of genes in a signature affect performance. We find that no matter how many genes were used in the signature, from 10–200, the Bayesian model consistently outperforms the frequentist.

(EPS)

Click here for additional data file.

Jaccard similarity of the top N genes, by statistical significance, between the Bayesian and frequentist meta-analysis methods, which is generally low.

(EPS)

Click here for additional data file.

Pathway analysis using Reactome PA—using the most significant genes for both methods– 1000 genes for Cardiomyopathy and Tuberculosis and 500 genes for Asthma.

A p-value cutoff of .2 was used for Cardiomyopathy, .2 for Asthma and .05 for Tuberculosis. BMA_only denotes genes only significant in Bayesian, but not frequentist. FMA_only denotes genes only significant in frequentist, but not Bayesian.

(EPS)

Click here for additional data file.

Datasets used for meta-analysis of asthma.

(XLSX)

Click here for additional data file.

Datasets used for meta-analysis of cardiomyopathy.

(XLSX)

Click here for additional data file.

Datasets used for meta-analysis of tuberculosis.

(XLSX)

Click here for additional data file.

27 Feb 2022

Dear Dr. Khatri,

Thank you very much for submitting your manuscript "Increasing reproducibility, robustness, and generalizability of biomarker selection from meta-analysis using Bayesian methodology" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Please address the concerns and comments raised by the reviewers; specifically, please clarify the choice of prior distributions and hyperparameters and assess their effects on the overall performance.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Yuchao Jiang, Ph.D.

Guest Editor

PLOS Computational Biology

Jian Ma

Deputy Editor

PLOS Computational Biology

***********************

Please address the concerns and comments raised by the reviewers; specifically, please clarify the choice of prior distributions and hyperparameters and assess their effects on the overall performance.

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Enclosed is a review of Kalesinskas et al’s manuscript

“Increasing reproducibility, robustness, and generalizability of biomarker selection from meta2 analysis using Bayesian methodology.”

In this paper, the authors describe a Bayesian meta-analysis approach that, in comparison to traditional frequentist approaches, is robust to outliers, increase the informativeness of between-study heterogeneity estimates, and reduces false positives and negatives. Overall, the paper is well-written and clearly motivated and explained. There are few comments I have that I detail section-by-section below. In summary, I think there may be room for added clarity for some figures.

INTRODUCTION

1. I commend the authors for clearly articulating this problem both from the statistical and biological perspectives.

2. In the last sentence of the Introduction, I suggest the authors name their method and R package.

RESULTS

1. The forest plots (like in Figure 2E) are hard to read and there is limited information about them in the figure captions. For example, what do the differently sized blue boxes mean?

2. I don’t quite see how effect sizes for FANCI are in opposite directions. The text seems to suggest the effect sizes are in opposite directions for the meta-analysis, but Figure 2E only shows different directions for one dataset.

3. The asterisks in Figure 4 are not defined. Are these P-value adjusted for multiple comparisons?

METHODS

1. How were the minimally informative prior hyperparameters elicited here? I.e., the Normal(0,3) and Uniform(0,2)?

Reviewer #2: In this study, the authors propose a Bayesian methodology for gene expression meta-analysis and show strong evidence that it outperforms standard frequentist approaches, and is, in particular, much less sensitive outlying data points and more sensitive to between-dataset heterogeneity. Overall, the study is well designed and the paper is well organized and written. I have no major concerns.

While I have no reservations about the final results, I think it would be helpful to include sensitivity analyses of the Bayesian prior parameters to assure that they are indeed non-informative, i.e. that the final gene rankings and posterior probabilities do not change by widening the prior further. I fully expect this to be the case.

Reviewer #3: The review is uploaded as an attachment.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: No: Authors provided a Github repo for the software (https://github.com/Khatri-Lab/bayesMetaIntegrator), but the codes and data to reproduce the analyses in this manuscript have not been provided.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: reviewer_comments.docx

Click here for additional data file.

9 May 2022

Submitted filename: Response to reviewers comments_V4.docx

Click here for additional data file.

29 May 2022

Dear Dr. Khatri,

We are pleased to inform you that your manuscript 'Increasing reproducibility, robustness, and generalizability of biomarker selection from meta-analysis using Bayesian methodology' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Yuchao Jiang, Ph.D.

Guest Editor

PLOS Computational Biology

Jian Ma

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: My comments have been addressed; I have no further comment. I congratulate the authors on their work and recommend publication!

Reviewer #2: The authors have addressed all of my comments.

Reviewer #3: The authors have answered all my questions. I don't have further questions.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: None

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Arjun Bhattacharya

Reviewer #2: No

Reviewer #3: No

22 Jun 2022

PCOMPBIOL-D-22-00035R1

Increasing reproducibility, robustness, and generalizability of biomarker selection from meta-analysis using Bayesian methodology

Dear Dr Khatri,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol