1.1 Motivation for evaluating biomarker sets

To motivate the need for evaluating association studies, consider first predictive studies, which use a labeled dataset, like the one shown at the top of Fig 1, to produce a predictive model (perhaps a decision tree, or a linear classifier) that can be used to classify future subjects–here, into one of the two classes: Y or N. In addition to the learned classifier, the researchers will also compute a meaningful estimate of its quality–i.e., of the accuracy (or AUROC, Kappa Score, etc.) of this classifier on an independent hold-out set [30], or the results of k-fold cross-validation over the training sample.

By contrast, many association studies report only a set of purported biomarkers, but provide no falsifiable claim about the accuracy of these biomarkers. Many meta-reviews claim that a set of biomarkers is problematic if they are not reproduced in subsequent studies [16, 31–33]. Given that biomarkers should be reproducible, we propose evaluating a biomarker set based on its reproducibility score. An accurate estimate of this score can help in at least the following three ways:

1. Researchers can compare various different “comparable” biomarker discovery algorithms to see which produces the biomarker set that is most reproducible. Here, “comparable” corresponds to the standard practice of only considering discovery tools that impose the same criterion, such as the same p-value, or only considering features that exhibit the same minimum fold-change. This type of analysis may help to determine errors within the biomarker discovery process.
   Moreover, we will see that MCC-correction, while useful in removing false-positives, can be detrimental to the goal of producing reproducible biomarkers; similarly, there is no reason to insist on p < 0.05 for the statistic test used.
2. A low reproducibility score suggests that few of the proposed biomarkers will be found in another dataset, and highlights the potential that these proposed biomarkers may not be accurate. This could motivate researchers to consider a dataset that is larger, to focus on a more homogenous population, or perhaps consider another biomarker discovery technique.
3. Finally, there are many meta-reviews [21, 34, 35] that note the lack of repeatability in many biomarker discovery papers, and question whether the techniques used are to blame. One way to address this concern is to require that each published paper include both the purported set of biomarkers, and also an estimate of its reproducibility score. The same way a prediction study’s “5-fold cross validation” accuracy tells the reader how accurate the classification model should be on new data, this biomarker-discovery reproducibility score will inform the reader whether to expect another study, on a similar dataset, will find many of the same biomarkers. Note that we should view the reproducibility score as necessary for considering a proposed model, but not sufficient–i.e., it might rule-out a proposed discovery model, but should not be enough to rule-in a model.

For these reasons, we provide an easy-to-use, publicly available webapp https://biomarker.shinyapps.io/BiomarkerReprod/ that anyone can use to produce meaningful estimates of the reproducibility of a set of biomarkers. (The underlying code is also available, from https://github.com/amirfrz/BMDA).

1.2 Related work

There have been many pairs of studies that have each produced biomarkers for the same disease or condition, but found little or no overlap between the two lists of purported biomarkers. Many papers have discussed this issue–some describing this problem in general [16, 33, 35], and others exploring specific examples [8, 32]. These papers suggest different reasons for the problem, such as the heterogeneous biological variations in some datasets [16, 33] or problems in the methods used that may lead to non-reproducible results [36, 37].

In particular, Zhang et al. [33] challenge the claim that the non-reproducibility problem in microarray studies is due to poor quality of microarray technology, by showing that inconsistencies occur even between technical replicates of the same dataset. They also show that heterogeneity in cancer pathology would further reduce reproducibility.

Ein-Dor et al. [16] also show the inconsistencies between the results of subsamples of a single dataset, demonstrating that the set of (gene) biomarkers discovered is not unique. They explain that there are many genes correlated with the group outcomes, but the empirical correlations change for different (sub)samples of instances. These two papers motivate our need for tools that can effectively estimate the reproducibility–such as the ones presented here.

Several projects [35–37] have attempted to formally analyse this problem. Ein-Dor et al. [35] describe a method, Probably Approximately Correct (PAC) sorting, that estimates the minimum number of instances needed for a desired level of reproducibility. As an example, this worst-case analysis proves that, to guarantee a 50% overlap between different gene lists for breast cancer, each dataset needs to include at least several thousand patients. This suggests poor repeatability results when using small sample sizes, which is consistent with our results for datasets with smaller sample sizes; see Subsection 3, especially Fig 4.

The goal of the MicroArray Quality Control (MAQC) project [36] was to address the problems and uncertainties about the microarray technology that were caused by the observation that different studies (of the same phenotype) often found very different biomarkers.

They suggest that the common approach of using just t-test p-values (especially with stringent p-values) can lead to poor reproducibility, which motivated them to consider methods like fold-change ranking with a non-stringent p cutoff, which they demonstrate leads to more reproducible gene sets. In a follow-up, Guo et al. [37] found similar results by using the same procedures for another dataset. However, Klebanov et al. [38] later show that these MAQC project results do not prove that using t-tests is necessarily unsuitable–i.e., just because another method (here fold-change) can generate more reproducible results, does not mean that it is performing better; as an extreme, the algorithm that declares every gene is a biomarker (think p = 1.0), is completely reproducible. They demonstrate these points by using a set of simulation studies (where they know the “true biomarkers”), and use either t-test or fold-change to propose potential biomarkers. These studies found that the t-test approach performed much better than the fold-change, in terms of recall (sensitivity). These results motivated us to use the t-test approach (rather than fold-change) as our main BD algorithm–which we use for all of our empirical experiments.

Our approximation algorithms use a type of re-sampling to bound reproducibility. Below we summarize several other studies that similarly deal with re-sampling and biomarker discovery. Some studies provide ways to better estimate the true statistical significance, but do not provide a framework for evaluating empirical reproducibility—e.g., Gagno et al. [39] used bootstrap resampling to estimate 95% confidence intervals and p-values for an internal assessment of their findings (related to breast cancer survival), Chitpin et al. [40] proposed a resampling-based method to better estimate the false discovery rate in chromatin immunoprecipitation experiments, and Pavelka et al. [41] proposed a resampling-based hypothesis testing algorithm that provides a control of the false positive rate for identification of differentially expressed genes. Furthermore, Alshawaqfeh et al. [42] and Zhao and Li [43] suggested methods for consistent biomarker detection in high-throughput datasets where evaluation was based on common biomarkers among the two resampled sets–i.e., they are considering the false positives but not false negatives. (By contrast, our use of Jaccard score involves both.) Other studies had different goals–e.g., Ma et al. [44] used resampling in a permutation test, to evaluate the predictive power of the identified gene set based on the accuracy of down-stream classification task; recall however that our goal is to evaluate the reproducibility of the biomarkers directly (not downstream). Filosi et al. [45] propose methods for evaluating the stability of reconstructed biological networks in terms of inference variability due to data subsampling; we however are focusing on the reproducibility of the individual biomarkers. Hua et al. [46] conducted a simulation study to compare the ranking performance of several gene set enrichment methods; by contrast, our approach considers the SET of biomarkers, not the ranking, and is over several real-world datasets (not just simulated ones). Note that none of these used re-sampling techniques to bound the expected replicability of the set of biomarkers found by some discovery algorithm, nor to demonstrate the validity of those bounds.

2.1 Formal description

As illustrated in Fig 1, a “Biomarker Discovery” algorithm, BD(⋅), takes as input a dataset D of n subjects, each described by r features F = {f₁, …, f_r} and labeled with a binary outcome, and returns a subset F′ ⊂ F of purported biomarkers.

Typically, each f ∈ F′ differs in some significant way between each class. To be more precise, let be the value of the i^th feature of the j^th subject, and ℓ^(j) be the outcome of the j^th subject (which is either + or -). Then a class difference means that the set of values of the i^th feature of the diseased individuals is significantly different from the values of that feature over the healthy individuals, . For simplicity, we will assume that these values are either all continuous (such as height, or the expression value of a gene), or all discrete (such as gender, or the genotype of a SNP). Subsection 2.2 below will describe several such biomarker discovery algorithms.

As noted in Eq 1, the Reproducibility Score RS(D, BD) quantifies the “reproducibility” of the set of proposed biomarkers BD(D) corresponding to running the biomarker discovery algorithm over the labeled dataset D. Here, we assume that the values of each feature , for each outcome c, are generated independently from a fixed distribution (i.e., “i.i.d.”)

Here and in general, we use P (⋅) to refer to either a probability density for continuous variables, or a probability mass for discrete variables. Note these are just the marginal distributions: we do not assume that the various features are independent from one another.–i.e., this does not necessarily correspond to Naive Bayes [30]. We will view as the matrix of these r × 2 different distributions, and let be the distribution for sampling n = n₊ + n₋ instances independently from this distribution, where instances are drawn from the distribution [p_1,+(⋅), …, p_r,+(⋅)] associated with positive outcomes, and similarly instances from the distribution [p_1,−(⋅), …, p_r,−(⋅)] associated with negative outcomes. Then for datasets sampled independently, we define (2) where the Jaccard score of two sets (3) is the ratio of the intersection to the union of these sets—hence J(A, B) ranges from 0 to 1, and is 1 if and only if A = B ≠ {}, and is 0 if and only if these sets are disjoint. (We define this to be 0 if A = B = {}.) Note that the Jaccard score is only one possible measure to evaluate the degree of overlap of gene signatures; Section C.4 in S1 Appendix discusses a slightly different measure that is sometimes used to measure the reproducibility of a set of biomarkers. See also Shi et al. [47] for a comprehensive overview of these measures.

Of course, we do not know , and so we use an empirical distribution , determined based on context from the dataset D, to produce the approximation (4) that estimates the reproducibility of the biomarker set BD(D), where the notation 〈D〉 = [|D⁺|, |D⁻|] refers to the pair of sizes of the positive and negative subjects in D—corresponding to [n₊, n₋]. Note that this reproducibility score deals with the sets of biomarkers that are produced by the BD(⋅) function, and not any specific biomarker.

Of course, Eq 4 suggests the obvious bootstrap sampling algorithm [48]. Empirically, however, we found that it did not perform well (see Section C.1 in S1 Appendix)–motivating the algorithms described in Subsection 2.3.

2.2 Biomarker discovery algorithms: BD(⋅)

We now discuss various approximation algorithms for biomarker discovery. As our goal is to illustrate the reproducibility issues with respect to the standard approach, we focus on that standard approach: where the biomarker discovery is based on independent two-sample t-tests, perhaps with some multiple comparison correction. This use of t-tests implicitly assumes that the feature values, for each outcome, are normally distributed. However, it is easy to adapt these algorithms to use other statistical tests, that do not make this distributional assumption.

See also Limitations in Section 4.

Recall that we are considering two types of datasets, depending on whether its feature values (the mentioned above) are continuous or discrete. However, for datasets with categorical values—SNPs in our analysis–we use a simple preprocessing step, which precedes all the BD(⋅) algorithms described here, to convert each categorical value to a real number: here converting each SNP feature, which ranges over the values { AA, Ab, bb }, to the real-values { 0, 1, 2 }, corresponding to the number of minor alleles (“b”) in the genotype. This allows us to view each such dataset as one with continuous values.

We assume that the real values of each feature, for each outcome, follows a normal distribution, which might be different for the different outcomes, and so we use an independent two-sample t-test for all of our empirical experiments. Recall that the test statistic is given by (5) where n₊ and n₋ are the number of instances with positive and negative outcomes, respectively, and with empirical means and and empirical variances and .

Note the biomarker discovery process essentially performs a single statistical test for each feature. As the number of features is often large–often tens-of-thousands, or more–many projects sought ways to reduce the chance of false discoveries; a standard way to do this is through some MCC method. We therefore consider biomarker discovery algorithms described as BD_{t,p, χ}(D), where the t in the subscript refers to the 2-sided t-test, the p for the p-value used, and χ to the MCC method. Our canonical example is BD_{t, 0.05, BH}(D), with p = 0.05, and χ = BH to the Benjamini/Hochberg correction [49]. This notation makes it easy to consider many variants –- e.g., adjusting the p-value used for the statistical test, whether it is applying another multiple testing correction, or none, etc. See Section C.2 in S1 Appendix for more details.

2.3 Algorithms that approximate the reproducibility score

As we have the dataset D with n ≈ [n₊, n₋] labeled instances, we can directly compute BD(D). To compute RS(D, BD(⋅)), however, we also need to produce one (or more) similar datasets D′, each with [n₊, n₋] subjects drawn from the same (implicit) distribution that generated D (with the same number of positive and negative instances), but which is presumably disjoint from D. While we do not have such D′’s, and so cannot directly compute the Reproducibility Score, we show below how to compute an overbound and an underbound of RS(D, BD(⋅)).

2.3.1 Overbound: oRS.

The oRS(D, BD(⋅), k) procedure produces (an estimate of) an overbound for RS(D, BD(⋅)), by making it easier for a feature to be selected to be in both purported biomarker sets. In this algorithm, D is our given fixed dataset, BD(⋅) is the given biomarker discovery algorithm, and k is a parameter to determine the number of trials when computing the overbound. (Here, and below, see the Glossary in Section A in S1 Appendix for a summary of the algorithms and their arguments.) The oRS algorithm first defines a size-2n dataset DD that contains two copies of each subject in D, of course with the same outcome both times. It then randomly partitions this DD into two disjoint size-n datasets D₁ and D₂, balanced by outcome. Here, each partition is with respect to the list of elements, so it will include duplicates. To insure that the resulting datasets are balanced, oRS first split DD into DD⁺ and DD⁻, where DD⁺ are the cases and DD⁻ the controls. It then forms by randomly drawing 1/2 of DD⁺, and by randomly drawing 1/2 of DD⁻, then merges ; see Fig 2. The dataset D₂ is then formed from the remaining subjects of DD not included in D₁. oRS then runs BD(⋅) on the datasets D₁ (resp., D₂) to produce two sets of biomarkers, and computes the Jaccard score for this pair of biomarker sets: J(BD(D₁), BD(D₂)). It then repeats this double-split-BD-Jaccard process k times, then returns the average of these k values: (6) where each dataset pair [D1_r, D2_r] is created independently using the above procedure.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Computing oRS.

Diagram showing oRS’s process of generating pairs of subsets for a dataset D (k times), then using those to compute oRS(D, BD(⋅), k).

https://doi.org/10.1371/journal.pone.0252697.g002

For each r, as we expect D1_r to overlap with D2_r, it is relatively likely that any D1_r-biomarker will also be a D2_r-biomarker (more likely than if D1_r was disjoint from D2_r), which means we expect the associated Jaccard score to be higher. This follows from the heuristic that, as two datasets have more common elements, we expect the number of biomarkers common to two datasets, to increase–i.e., (7)ceteris paribus. Here, we view each of {A₁, A₂, B₁, B₂} as a set of n = n₊ + n₋r-dimensional instances. Note this relationship is simply a heuristic to motivate the algorithm–one that we expect to hold in practice. Section D in S1 Appendix provides some basic arguments, and empirical evidence, to support this claim.

We close with three quick observations:

1. Expected Overlap: We expect 50% of the instances to be duplicated in any given pair of datasets. See Lemma 2 in Section D.4 in S1 Appendix.
2. Relation to Bootstrap Samples: Here, each subject occurs exactly twice in each pair of datasets D1_r and D2_,r. If we instead used bootstrap sampling (called bRS below), we expect many subjects would occur more often in the pair of datasets. (See Lemma 3 in Section D.4 in S1 Appendix.) Given Heuristic 7, this means we expect bRS’s Jaccard score here to be higher than for oRS’s; as oRS is already an overbound for the true score (RS), this means bRS would be a worse bound, as bRS ≥ oRS ≥ RS. This is why we use our “doubling approach” oRS rather than bootstrap sampling bRS, as oRS produces values that are smaller, but still remains an overbound, as desired. See Figs 7 and 8 and Section C.1 in S1 Appendix.
3. Relation to RS*: Our oRS approach is clearly related to RS* (Eq 2), as both compute the average Jaccard score of pairs of size-n datasets sampled from a distribution. They differ as (a) RS* uses the true distribution , while oRS uses only the estimate , (b) RS* is the true average while oRS is just the empirical average over k trials, and (c) RS* will draw independent datasets, but the oRS datasets will overlap.

2.3.2 Underbound: uRS.

The uRS(D, BD, k) procedure produces (an estimate of) an underbound for RS(D, BD(⋅)) by making it harder for a feature to be selected to be in both purported biomarker sets. First, observe that as [n₊, n₋] increases (keeping the n₊-to-n₋ ratio fixed, as we consider changing the size of the dataset), we expect the statistical estimates to be more accurate, and in particular, statistical tests for differences between the two classes will be correct more often. Hence, a statistical test will better identify the “true” biomarkers F* from a size-n subset D⁽ⁿ⁾, versus from a size-n/2 subset D^(n/2). Now consider two size-n datasets and , and also two size-n/2 datasets and . As and are each closer to F* than and , we expect and to be closer to each other, than and , which means we expect that

In general, given that (8) (where is the empirical average over k samples), this argues that the RS* score should increase with the size s of the dataset–which suggests that

Fig 4(a) presents empirical evidence, over 5 datasets, supporting this claim —showing that the Jaccard score increases as we increase the size s of the datasets. Section D in S1 Appendix provides some arguments, and additional empirical evidence (over hundreds of simulations), that further support this heuristic.

This motivates our underbound algorithm uRS(D, BD, k), which first partitions D into two disjoint size-n/2 subsets, E₁ and E₂, with balanced outcomes. It then computes J(BD(E₁), BD(E₂)) which, assuming Heuristic 9, is an underbound in expectation for RS(D, BD(⋅)). uRS does this partitioning k times, producing k different dataset pairs {[E_1,r, E_2,r]}_{r = 1, …, k}, and returning the average Jaccard score, i.e. (10)

See Fig 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Computing uRS.

Diagram showing uRS’s process of generating pairs of subsets for a dataset D (k times), then using those to compute uRS(D, BD(⋅), k).

https://doi.org/10.1371/journal.pone.0252697.g003

2.4 Empirical study over various datasets

There are now many publicly-available datasets that have been used in association studies. Here, we use them to …

1. U1: Better understand what Jaccard scores are typical, for the standard BD(⋅) algorithms;
2. U2: Determine whether our predictions match the results of earlier meta-analyses; and
3. U3: Determine if our approximations are meaningful—i.e., if (for reasonable values of k): (11)(12)

The next section will explicitly discuss all three issues. Of course, given only a single dataset D of size-n, we cannot compute, nor even estimate, the true value of RS(D, BD(⋅)). However, we can estimate RS(D^(n/2), BD(⋅)), where D^(n/2) is a size-n/2 (outcome-balanced) subset of D. In fact, uRS(D, BD(⋅), k) is a meaningful estimate of RS(D^(n/2), BD(⋅)); below we will use (13)

We will then compare this against uRS(D^(n/2), BD, k) and oRS(D^(n/2), BD, k) and to see whether the relations of Eqs 11 and 12 both hold, with respect to various size-n/2 subsets D^(n/2).

More generally, we can do this for any size-s subset D^(s) of D where s ≤ n/2. Here, we need a set of pairs of disjoint outcome-balanced subsets D′, D′′ ⊂ D where |D′| = |D′′| = s and D′∩D′′ = {}. For a fixed dataset D, and specified number , we can then plot these values along with oRS(D^(s), BD, k) and uRS(D^(s), BD, k), as a function of s to see their behaviour; see Fig 4(b), for the Metabric dataset. Our website https://biomarker.shinyapps.io/BiomarkerReprod/ also provides this visualization.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Reproducibility Scores for different subset sizes.

(a) For each of the 5 datasets D, each point shows the average (± sd) Jaccard over k = 20 pairs [D1^(s), D2^(s)] of disjoint size-s subsets of D. Here, n is the size of the original dataset–note this can only go to n/2–and we are using the standard BD_{t, 0.05, BH}. (b) Showing how the approximations relate to one another, and scale with the size s of the dataset. Here we are using subsets of the Metabric dataset, with n = 1654. We observed the same behavior for all datasets.

https://doi.org/10.1371/journal.pone.0252697.g004

We explored our approximations over 25 different real-world datasets, including 16 microarray datasets and 2 RNAseq datasets with continuous data, and 7 SNP datasets with categorical data (see Table 1). This first set includes 4 of the gene expression datasets discussed in the Zou et al. [8] meta-analysis—each describing metastatic versus non-metastatic breast primary cancer subjects—to see if our method is consistent with their empirical results. We also included 11 other relatively-small gene expression datasets (with 19 to 187 subjects), focusing on human studies that had a binary class outcome from the GEO repository. To explore how our tools scale with size, we also included 3 other relatively large datasets, with 532 to 1654 subjects. As these were survival datasets, we set the binary outcome based on the median survival time (removing any subject that was censored before that median time). In addition to these 4+11+3 = 18 gene expression datasets (with real-valued entries), we also include 7 SNP datasets (with 39 to 164 subjects), with discrete values, also selected from human studies with binary class outcome s. Fig 5 plots the number of features and biomarkers found, using the BD_{t, 0.05, BH} algorithm, for each dataset—both D⁽ⁿ⁾ and D^(n/2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Results for all 25 datasets when using all the subjects.

This table is sorted by sample size (#subjects)–corresponding to Fig 5. The first 18 entries are Gene Expression datasets (including the 4 “*”ed entries, from the Zou et al. [8] meta-study), and the final 7 are SNP datasets. Reproducibility Scores are shown in the form of mean ± standard deviation. The “(Majority %)” values are the percentage of the subjects in the dataset with the more common outcome–e.g., 53% of the subjects in the GDS968 dataset are labeled “+” (for “long survival time”), and 82% of GSE7390 are labeled “−” (for “Non-Metastatic”).

https://doi.org/10.1371/journal.pone.0252697.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Number of biomarkers found for D^(n/2) compared to D⁽ⁿ⁾.

Box+whiskers plots showing number of biomarkers found for D^(n/2) when using BD_{t, 0.05, BH} over k = 20 iterations for various dataset, compared to the number of biomarkers for D⁽ⁿ⁾, and to the number of features in each dataset. Note the y-axis is a log-scale. (Note we first changed all “0” values to “10⁻¹”.) For details, see Tables 1 and 2.

https://doi.org/10.1371/journal.pone.0252697.g005

4.1 Recommended use

As suggested above, whenever researchers have identified a possible set of biomarkers from a dataset D and Biomarker Discovery Algorithm BD, we encourage them to apply our analytic technique, and where appropriate, even our specific bounding algorithms (from https://biomarker.shinyapps.io/BiomarkerReprod/), to estimate the Reproducibility Score of this proposed set. If those scores (especially the lower bound uRS(D, BD)) are sufficiently high, they can use those biomarkers, confident that they (as a set) are reproducible. But if not, the researchers may want to explore other ways to identify reproducible biomarkers. One obvious approach is to use a larger dataset, with more instances, as we know that RS increases with the size of the dataset |D|; see Heuristic 7 and the analysis in Section D.2 in S1 Appendix. Alternatively (or in addition), recall the reproducibility depends on both the dataset, and the Biomarker Discovery Algorithm BD; perhaps some other BD_{test,p − val, MCC} would work better here? We could consider modifying all 3 of the components:

• While the MCC methods are designed to reduce the false positives, it is not clear whether they improve RS. Indeed, our experiments (Fig C.1 in S1 Appendix) show that RS values reduce with MCC—which points to using BD_{test,p, None}.
• The non-reproducibility problem might be due to the statistical test used. As noted earlier, the t-test implicitly assumes a Gaussian distribution of the features (or the normality of residuals); perhaps another test would work better for some dataset / distribution of instances.
• Finally, they may want to modify the p-value, p, as other values of p may lead to better reproducibility.

4.2 Limitations

This paper provides an effective way to estimate the reproducibility of the biomarkers found from a dataset, using a biomarker discovery algorithm. While the message of this paper is very general, the specific analyses here all used the standard discovery algorithm BD_{t, 0.05, BH}, to illustrate that the issues apply to the approaches commonly used. Section C in S1 Appendix explores some other discovery methods, that are also based on t-tests. While this use of t-tests implicitly assumes normality of the feature values, nothing in our general approach relies on this specific test—we could use other tests that make other parametric assumptions. (Note that our analysis appears to work effectively even when dealing with some Bernoulli (non-normal) features.) We anticipate the same approach would hold for other tests, such as Mann-Whitney, Wilcoxon or even multivariate analysis–but this is future work. All of our empirical studies dealt with standard datasets, whose outcomes are binary and whose values were either all real values or all categorical values; none had some of each. Our analytic model considers the overlap of biomarkers found from two datasets, of the same size—e.g., we do not consider how the biomarkers obtained from a size-100 dataset, overlap with those from a size-300 dataset. Finally, the recommendations of the previous subsection suggest another future direction: produce the “BD′(D, s) algorithm” that, given a dataset D and a minimum score s > 0, returns the parameters [test, p-val, MCC] for a biomarker discovery algorithm BD_{test,p − val, MCC} that would produce a biomarker set whose Reproducibility Score would be at least s–i.e., we expect that RS(D, BD_{test,p − val, MCC}(⋅)) ≥ s, suggesting this level of reproducibility for the proposed biomarkers BD_{test,p − val, MCC}(D).

4.3 Contributions

This paper has several contributions: (1) It motivates, then provides, a formal definition of reproducibility, that can help researchers evaluate a set of purported biomarkers; (2) It provides a pair of algorithms that can accurately bound this “reproducibility score” for a given dataset and biomarker discovery algorithm; (3) It provides empirical evaluation of these algorithms, over 25 different real-world datasets, to demonstrate that they work effectively; and (4) It introduces a freely-available website https://biomarker.shinyapps.io/BiomarkerReprod/ that runs these algorithms on the dataset entered by a user, and biomarker discovery system, which will allow users to easily evaluate the quality of the biomarker set produced. Given how easy it is to use this tool, we hope that future researchers will automatically use it to quickly evaluate the quality of the purported biomarkers, then include these estimates when they publish their biomarkers. We anticipate this analysis may also lead to new biomarker discovery algorithms, to optimize this reproducibility score.