Author response: Extent, impact, and mitigation of batch effects in tumor biomarker studies using tissue microarrays (2021)

Article Figures and data Abstract Editor's evaluation eLife digest Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Tissue microarrays (TMAs) have been used in thousands of cancer biomarker studies. To what extent batch effects, measurement error in biomarker levels between slides, affects TMA-based studies has not been assessed systematically. We evaluated 20 protein biomarkers on 14 TMAs with prospectively collected tumor tissue from 1448 primary prostate cancers. In half of the biomarkers, more than 10% of biomarker variance was attributable to between-TMA differences (range, 1–48%). We implemented different methods to mitigate batch effects (R package batchtma), tested in plasmode simulation. Biomarker levels were more similar between mitigation approaches compared to uncorrected values. For some biomarkers, associations with clinical features changed substantially after addressing batch effects. Batch effects and resulting bias are not an error of an individual study but an inherent feature of TMA-based protein biomarker studies. They always need to be considered during study design and addressed analytically in studies using more than one TMA. Editor's evaluation Tissue microarrays (TMA) have become a mainstay in clinical and basic research, for both discovery and validation of biomarkers. This manuscript provides relevant methodological considerations for cancer researchers investigating tissue-biomarkers using TMAs. A comprehensive investigation was conducted using a combination of analytic approaches using empirical data and simulated data to support key findings and conclusions. The authors approach the possible sampling variation in a thoughtful way, not only quantifying the issue systematically, but working towards a solution. https://doi.org/10.7554/eLife.71265.sa0 Decision letter Reviews on Sciety eLife's review process eLife digest To understand cancer, researchers need to know which molecules tumor cells use. These so-called 'biomarkers' tag cancer cells as being different from healthy cells, and can be used to predict how aggressive a tumor may be, or how well it might respond to treatment. A popular technique for assessing biomarkers across multiple tumors is to use tissue microarrays. This involves taking samples from different tumors and embedding them in a block of wax, which is then cut into micro-thin slices and stained with reagents that can detect specific biomarkers, such as proteins. Each block contains hundreds of samples, which all experience the same conditions. So, any patterns detected in the staining are likely to represent real variations in the biomarkers present. Many cancer studies, however, often compare samples from multiple tissue microarrays, which may increase the risk of technical artifacts: for example, staining may look stronger in one batch of tissue samples than another, even though the amount of biomarker present in these different arrays is roughly the same. These 'batch effects' could potentially bias the results of the experiment and lead to the identification of misleading patterns. To evaluate how batch effects impact tissue microarray studies, Stopsack et al. examined 14 wax blocks which contained tumor samples from 1,448 men with prostate cancer. This revealed that for some biomarkers, but not others, there were noticeable differences between tissue microarrays that were clearly the result of batch effects. Stopsack et al. then tested six different ways of fixing these discrepancies using statistical methods. All six approaches were successful, even if the arrays included tumors with different characteristics, such as tumors that had been diagnosed more or less recently. This work highlights the importance of considering batch effects when using tissue microarrays to study cancer. Stopsack et al. have used their statistical approaches to develop freely available software which can reduce the biases that sometimes arise from these technical artifacts. This could help researchers avoid misleading patterns in their data and make it easier to detect real variations in the biomarkers present between tumor samples. Introduction Tissue microarrays (TMAs) were first developed in the 1990s as an efficient way to examine tissue-based biomarkers (Kononen et al., 1998). Since then, TMAs have been used in thousands of studies to evaluate histologic and molecular biomarkers, mostly in cancer tissue. Individual TMAs consist of cylindrical cores from hundreds of tissue samples embedded in one recipient block (Kononen et al., 1998; Kallioniemi et al., 2001). Studies often include more than one TMA. Even when biomarker assays are well standardized and run conditions are diligently kept fixed, some TMA slides (batches) may have measurements systematically too low or too high, and some batches may have wider spread around the true values of the biomarker than others. In general, such batch effects can have a profound impact on the validity of biomarker studies, such as those using RNA microarrays (Tworoger and Hankinson, 2006; Leek et al., 2010). Contrary to popular belief, whether such measurement error induces upward or downward bias in results is not guaranteed to follow simple heuristics (van Smeden et al., 2020). Whether and to what extent TMAs are affected by batch effects has not been empirically assessed. TMAs pose unique challenges. For example, when tumor tissue is collected prospectively for inclusion on TMAs, tumor characteristics may differ between batches due to nonrandom assignment of cases, as well as temporal trends in tumor risk factors, screening, and diagnosis. Differences in tissue processing or storage across tissue specimens may have differential impact on biomarkers. Including calibration samples for quality control is also more challenging for TMAs than, for example, assaying of blood samples, because repeat sections from a tumor may differ due to intratumoral heterogeneity rather than only batch effects. In this study, we assess batch effects in a large set of centrally constructed TMAs from prostate cancer tissue from 1448 men in two nationwide cohort studies. We quantify the extent to which protein biomarker variation could be explained by batch effects. We probe different methods for mitigating batch effects while maintaining true, "biological," between-TMA variation, including in a plasmode simulation. Finally, we demonstrate the impact of handling batch effects on commonly performed biomarker analyses. Results Extent and type of batch effects To evaluate the presence of batch effects in studies using TMAs, we studied tumor tissue from 1448 men with primary prostate cancer on 14 TMAs (labeled "A" through "N"), each including multiple tumor cores from 47 to 158 patients per TMA (Figure 1). Multiple cores from the same tumor (usually 3) were always located on the same TMA. Figure 1 Download asset Open asset Characteristics of men with prostate cancer with tissue included on the 14 tumor tissue microarrays. (A) Calendar years of cancer diagnosis, with thick lines indicating median, boxes interquartile ranges, and whiskers 1.5 times the interquartile range. (B) Counts of tumors by Gleason score. (C) Counts of tumors by pathological TNM stage (RP: radical prostatectomy). (D) Rates of lethal disease (metastases or prostate cancer-specific death over long-term follow-up), with bars indicating 95% confidence intervals. As throughout, multiple cores are summarized per tumor. TMAs were used to quantify 20 protein biomarkers (Figure 2). Biomarker values showed noticeable between-TMA variation, despite immunohistochemical staining having been conducted at the same time for all 14 TMAs. We estimated that across the 20 biomarkers, between-TMA variation explained between 1% and 48% of the overall variation in biomarker levels (intraclass correlation coefficient, ICC), with half of the biomarkers having ICCs greater than 10% (Figure 2). Figure 2 with 1 supplement see all Download asset Open asset Biomarkers stained, staining and scoring methods, and intraclass correlation coefficients (ICCs). (A) Tissue microarrays (TMAs) assessed for each marker (dark blue, yes). (B) Approach to staining biomarkers: automated staining system versus manual staining (gray, yes); quantification of biomarkers: software-based scoring versus eye scoring (blue, yes); biomarker quality assessed: staining intensity, proportion of cells positive for the biomarker, area of tissue positive for the biomarker (yellow, yes). (C) Counts of tumors assessed for each biomarker. (D) Between-TMA ICCs (i.e., proportion of variance explained by between-TMA differences) for each biomarker, with 95% confidence intervals. Empty symbols indicate the 97.5th percentile of the null distribution of the ICC obtained by permuting tumor assignments to TMAs; asterisks indicate between-Gleason grade group ICCs. Biomarkers are arranged by descending between-TMA ICC. In an example biomarker, estrogen receptor alpha in nuclei of stromal cells (Figure 3), the means of the most extreme TMAs differed by 2.2 standard deviations in intensity of expression and variances differed up to 9.3-fold. Other biomarkers showed similar between-TMA variation by magnitude and by which TMAs had the most extreme values (Figure 4A). Likewise, we observed that not only means, but also variances of biomarker levels differed between TMAs, although patterns of heteroskedasticity appeared weaker than for means (Figure 4—figure supplement 1). In contrast, we found little evidence for more complex patterns of batch effects, such that tumors with specific grade, stage, or year of diagnosis would have been particularly affected by between-TMA differences (Supplementary file 1a). Nevertheless, observations from the same TMAs tended to be clustered together when projected onto the first two principal components, capturing 27% of the variance in all biomarkers (Figure 4B). Figure 3 with 1 supplement see all Download asset Open asset Effect of batch effect mitigation on a biomarker with strong between-tissue microarray (TMA) variation. (A) The protein biomarker estrogen receptor-alpha was quantified as staining intensity in nuclei of epithelial cells, averaged over all cores of each tumor. Each panel shows processed data for a specific approach to correcting batch effects. Notes in the upper right corner indicate which properties of batch effects were potentially addressed. Each data point represents one tumor. y-axes are standard deviations of the combined data for the specific method. Thick lines indicate medians, boxes interquartile ranges, and whisker length is 1.5 times the interquartile range. (B) Example photographs of TMAs; brown color indicates positive staining. Figure 4 with 7 supplements see all Download asset Open asset Patterns, source, and remediation of batch effects. (A), Biomarker mean levels by tissue microarray (TMA), in biomarker-specific standard deviations (y-axis). Each point is one TMA. (B) First two principal components of biomarkers levels on all 14 TMAs, with color/shape denoting TMA. Each point is one tumor. (C) The same first two principal components, with color/shape denoting Gleason score. (D) Per-core biomarker levels for tumors with multiple cores included on two separate TMAs, for estrogen receptor (ER) alpha and beta, both in standard deviations. Each point is one tumor core. (E) Pearson correlation coefficients r between uncorrected and corrected biomarker levels. Entries are averages across all markers. Some biomarkers were stained using automated staining systems, other stains were done manually (Figure 2). Moreover, the method of scoring, including human (eye) scoring and computer-assisted quantification, differed between biomarkers, as did the main quantitative score, typically a measure of staining intensity, a proportion of cells above an intensity threshold, or a combination of both (Figure 2). Notably, between-TMA differences were present with any of these approaches. For example, batch effects were not only present when considering intensities of biomarker staining, as for the estrogen receptor alpha and beta example. Even when setting cutoffs for staining visible by eye and quantifying the number of stain-positive cells, 8% (95% confidence interval [CI], 2–15) of variance in estrogen receptor alpha positivity and 27% (95% CI, 11–42) of estrogen receptor beta positivity were attributable to between-TMA variation (Figure 2—figure supplement 1). Our data do not allow distinguishing which of these approaches, if any, were less prone to batch effects. In summary, we observed a large and concerning degree of between-TMA variation for several biomarkers that were quantified using different approaches, suggesting that addressing batch effects could significantly impact scientific inference. Source of batch effects The noticeable proportion of variance attributable to TMAs could have two possibly co-existing explanations. First, that between-TMA differences in biomarkers reflect different patient and tumor characteristics that need to be retained. Second, that between-TMA differences are artifacts due to systematic measurement error that need to be removed (batch effects). In support of the first hypothesis, there were noticeable differences in patient and tumor characteristics between TMAs that are likely associated with biomarker levels (Figure 1). Along with a 14-year range between the per-TMA medians of cancer diagnosis year, there were differences in the proportion of tumors with a Gleason score of 8 or higher (between 11% and 33%) and rates of lethal disease (between 2 and 16 events per 1000 person-years of follow-up). In support of the second hypothesis, we observed that certain TMAs had consistently higher or lower biomarker values for the majority of tested biomarkers (Figure 4A). For example, the same batches that showed higher-than-average biomarker values for stathmin also had higher-than-average values for PTEN. This example is noteworthy because both markers were assayed together on the same section of each TMA using multiplex immunofluorescence, and stathmin would be expected to be expressed in more aggressive tumors with activation of the PI3K signaling pathway while PTEN expression would be expected to be low in the same tumors (Stopsack et al., 2020). Further supporting the second hypothesis, we did not observe any meaningful reduction in ICCs when we considered tumors that had the same clinical features in terms of Gleason score and stage (Figure 4—figure supplement 2). Moreover, the association between Gleason score and biomarker levels (Figure 2D) was considerably lower than between TMAs and biomarker levels, as underscored by less pronounced visual separation of principal components by Gleason score (Figure 4C) than by TMA (Figure 4B). Gleason score differences explained no more than 13% of the variance in biomarker levels (for prostate-specific membrane antigen, PSMA; 95% CI for ICC, 0.02–0.29), and 13 of the 20 biomarkers had ICCs by Gleason score of 1% or less (Figure 4—figure supplement 3). To directly disentangle both hypotheses, we further examined data on 10 tumors with a total of 53 tumor cores for which some cores were included on different TMAs (Figure 4D). These were not included in the previous analyses and had estrogen receptor scoring data. This design allowed us to estimate biomarker differences directly attributable to between-TMA variability within the same tumors while controlling for the between-core variability expected due to intratumoral heterogeneity. Of the total variance in estrogen receptor alpha levels, 30% (95% CI, 0–67) was explained by between-TMA variation; for estrogen receptor beta, 24% (95% CI, 0–60) was explained by between-TMA variation. For comparison, between-tumor variation explained 37% (95% CI, 4–68) of the variance of estrogen receptor alpha levels and 26% (95% CI, 0–57) of the variance of estrogen receptor beta levels. Collectively, while these observations highlighted moderate differences in clinical and pathological characteristics between TMAs, they suggested that between-TMA differences were largely due to batch effects. Mitigation of batch effects We implemented six different approaches for batch effects mitigation and compared these to the uncorrected biomarker levels (Figure 3, Figure 3—figure supplement 1). Two mitigation approaches, batch means (approach 2) and quantile normalization (approach 6), assumed no true difference between TMAs is arising from patient and tumor characteristics, while all other approaches attempted to retain such differences between TMAs. It is possible that the choice of mitigation approaches may be optimized using knowledge of the source of the batch effect. This would be the case if each method "specialized" in mitigating effect from specific sources. We have not investigated this possibility here. Overall, correlations between values adjusted by different approaches were higher (mean Pearson r, 0.97–1.00) than between uncorrected values and corrected values (r, 0.90–0.95), regardless of mitigation approach (Figure 4E). Approaches 2–7 reduced visible separation by batch on plots of the first two principal components (Figure 4—figure supplement 4). Variance attributable to between-TMA differences decreased to ICCs of <1% for all markers (Supplementary file 1b). An exception was the quantile regression-based approach 5; the ICCs after this approach remained up to 10%. This method does not explicitly address differences in means between batches but allows associations between clinical and pathological factors and biomarker levels to differ at high and low quantiles (Figure 4—figure supplement 5). The differences between uncorrected values and batch effect-corrected values were remarkably similar between the mean-based approaches using approaches 2 (simple means), 3 (standardized batch means), and 4 (inverse probability-weighted batch means; Figure 4—figure supplement 6). Consequently, batch effect-corrected values by approaches 2–4 were highly correlated (Figure 4E). All mean-only batch effect mitigations also gave the same results when fitting outcome models stratified by batch (Figure 4—figure supplement 7). However, batch-specific results differed for approaches that targeted between-batch differences in the variance of biomarkers. Validating batch effect mitigation in plasmode simulation To compare the performance of the different batch mitigation approaches in a time-to-event analysis, we applied plasmode simulation (Franklin et al., 2014) to fix the expected strength of the biomarker exposure–outcome relationship a priori before artificially introducing batch effects. The correlation structure between biomarker and confounders and between confounders and batches from the actual data (Figure 5—figure supplement 1A, C) was preserved in the plasmode-simulated data. Likewise, across a range of hazard ratios for the biomarker–outcome association, confounder–outcome associations remained unchanged (Figure 5—figure supplement 1B, D). We first evaluated a setting in which we did not introduce batch effects (Figure 5A). Here, the observed hazard ratios without batch effect mitigation equaled the expected. When performing (unnecessary) batch effect mitigation, observed hazard ratios were still comparable with the expected hazard ratios (Figure 5D; see Supplementary file 1c for CIs). Figure 5 with 3 supplements see all Download asset Open asset Plasmode simulation results. (A–C) Biomarker levels by tissue microarray in three simulation scenarios. (D–F) True versus observed hazard ratios for the biomarker–outcome association after alternative approaches to batch effect correction, with correction methods being numbered as in the Materials and methods section. The solid blue line indicates no correction for batch effects. (A, D) no batch effects; (B, E) means-only batch effects; (C, F) means and variance batch effects. We then introduced batch effects by adding batch-specific mean differences to the observed biomarker levels, yet without introducing differences in variance by batch (Figure 5B). Without batch effect mitigation, for a true hazard ratio of 3.0, the observed hazard ratio, averaged over simulations, was 2.17 (95% CI, 1.86–2.53), an underestimate by 28% (Figure 5E; Supplementary file 1c). In contrast, all mitigation approaches produced CIs that covered the expected hazard ratio (e.g., approach 6 quantile normalization: hazard ratio, 3.03; 95% CI, 2.48–3.69). When we introduced batch-specific differences in both means and in variances (Figure 5C), the observed hazard ratio without batch effect mitigation decreased to 1.90 (95% CI, 1.66–2.16) compared to the expected hazard ratio of 3.0 (Figure 5F; Supplementary file 1c). Batch effect mitigation methods that only focus on means (approaches 2–4) reduced but did not fully eliminate bias, with hazard ratios ranging between 2.67 and 2.70. Methods that address differences in both mean and variance resulted in less bias, with an observed hazard ratio of 3.11 (95% CI, 2.54–3.81) for approach 6 (quantile We also included two approaches. models by by (approach resulted in but was less efficient than other approaches, with a risk of bias, and resulted in considerably wider CIs in simulation. Including batch as a in a (approach was and A of both approaches is that they do not explicitly estimate batch biomarker values that could be evaluated were on the of confounders between batches and at most associations between the biomarker and resulting in We introduced both and strong associations between biomarker and confounders and more between batches (Figure 5—figure supplement 2). In all the of mitigation methods was preserved (Figure 5—figure supplement 3, Supplementary file with the bias obtained through quantile normalization (approach 6). when using uncorrected biomarker levels in the presence of any batch effects, if there was no association between biomarker and outcome (i.e., a hazard ratio of and with mean-only approaches 2–4 if variance was also affected by batch effects. In no possibly with the quantile regression-based approach were after batch effect mitigation from the expected values than results on uncorrected biomarker levels. of batch effects To how batch effect mitigations the results of commonly conducted tumor biomarker we estimated how uncorrected and corrected biomarker levels were associated with Gleason score and with rates of lethal For markers with little between-TMA variability such as there were no noticeable differences in associations between using and adjusted biomarker levels of as expected from plasmode simulation. However, for markers with higher between-TMA variability and stronger associations with the approaches to noticeable differences (Figure 6). For example, uncorrected stathmin expression levels were not associated Gleason score standard deviations per one 95% CI, to while the difference in levels corrected to approach 6 (quantile was (95% CI, to suggesting a potentially different (Figure Supplementary file In models for lethal disease (Figure the hazard ratio for the of the compared to the was (95% CI, after mitigation using approach the hazard ratio was (95% CI, suggesting that biomarker levels could underestimate the association by (Supplementary file Figure 6 Download asset Open asset of batch effect mitigation on scientific inference. (A) Gleason score and biomarker levels standard per Gleason grade (B) Biomarker levels and to lethal with hazard ratios per one standard increase in biomarker levels from In both blue indicate using uncorrected biomarker levels, indicate batch effect-corrected levels, quantile normalization (approach 6). are 95% confidence intervals. Biomarkers are by between-tissue microarray intraclass correlation Discussion The key strength of using TMAs is their in the of biomarkers on a large number of tissue specimens (Kononen et al., 1998). to other batch effects have to be considered in TMA biomarker As we for some of the biomarkers, batch effects can be of We that batch effect mitigation is possible and can study In study of prostate tumor between-TMA differences explained 10% or more of the variance in biomarker levels for half of the included biomarkers, considerably more than one of the pathological features in prostate cancer, Gleason All mitigation approaches to reduce batch effects, whether they attempted to retain real differences between tumors from different TMAs or to corrected biomarker levels that were more similar to each other than they in general, to the uncorrected biomarker levels. In from a large set of protein tumors biomarkers in prostate cancer, we how mitigating batch effects results and their validity for biomarkers affected by batch effects. batch effects between TMAs are when a how tumor samples are and assayed is as are stratified or of tumor samples to different TMAs (Tworoger and Hankinson, when However, the batch effects that we observed despite all Moreover, samples be collected and TMAs may be constructed in studies over were differences in the clinical and pathological characteristics between TMAs, an issue that may be in studies. of tumors to TMAs that for example by of lethal tumors with of may not be Likewise, in of the thousands of studies using TMAs it be possible to tumors to different TMAs and repeat all work to reduce the of batch effects. An in the design is that tissue samples are and be blood samples. tumor samples that could as a quantitative calibration for all biomarkers do not is to include lines that have been and on each TMA. lines address of the lines are often unique and as such may not be relevant for all biomarkers. approach is to include samples from the same tumor case across TMAs, which would allow for of batch effects. For these a approach that batch effects and them analytically is to batch effects during the study design we the when TMA-based tissue biomarker studies (Figure 7). First, the extent of batch effects be and in any study of cancer tissue using TMAs. TMA slides and plots (Figure et al., is Between-TMA variation be for example by that to variation of biomarker levels between TMAs compared to that between or within tumors and 2010). In study, for half of the biomarkers, ICCs for between-TMA variation were at less than although the proportion of batch variation be on the Whether TMAs differ in terms of biomarker levels, low levels of or variability between tumors also what impact of between-TMA differences to Figure 7 Download asset Open asset for and handling batch effects between tissue microarrays approaches at the design of a all TMA-based studies the for batch effects a biomarker has been batch effects only be there is no for their Batch effect-corrected biomarker levels can be by the Second, the source