BACKGROUND
there is a wealth of information enclosed in the massive amount of microarray data so far accumulated in public repositories.
the variety of data sets generated from the assortment of experiments is a major obstacle in the path leading from these data to information.
specific issues relating to data heterogeneity across microarray studies include differences across platforms, differences within and between labs, and/or differences among experimental factors such as treatments and tissues  <cit> .
furthermore, concerns regarding integration of studies from multiple sources in general, such as variations in design, research goals, or quality of implementation, add to these issues  <cit> .
inappropriate integration of microarray data from public repositories could lead to equivocal results  <cit> .
the “simpson’s paradox”  <cit> , which refers to contradictory statistical results obtained when analysis is performed within versus across groups of data  <cit> , is an example of mishandling of data.
blyth  <cit>  gives an example involving the analysis of 2x <dig> contingency tables across two groups, and hassler and thadewald  <cit>  also illustrate simpson’s paradox when correlation coefficients are estimated from a pool of two groups versus within each group.
in both cases, the paradox can be explained as results are further investigated in light of the specific statistical properties of each group of data.
the “ecological fallacy” happens when the correlation of aggregated variables results in a significant relationship that is due only to aggregation rather than to any real association  <cit>  .
an early example of an ecological fallacy can be found in gehlke and biehl  <cit> , whose study of grouping effects in census tract data showed that the magnitude of correlation coefficients of two variables tend to increase as the level of census tract aggregation increases.
this problem was later referred to as the “modifiable areal unit problem” and further studied by openshaw and taylor  <cit> .
combining statistical results  of independent studies that address similar questions has been a standard procedure in classic meta-analyses  <cit> .
this approach entails analyzing each data set independently and then combining the results, as in a mosaic.
meta-analysis of microarray data has been applied in a broader context, as some works include data spanning a wide range of purposes and designs.
parmigiani et al.
<cit> , in a quest for a common gene signature across multiple cancer types, developed a statistical method to identify and assess the intersection of multiple gene expression signatures across  <dig> published cancer-related microarray studies.
on the other hand, wirapati et al.
<cit>  and rhodes et al.
<cit>  developed specific meta-analysis methods to integrate gene expression signatures of breast and lung cancer, respectively, across independent studies of microarray data.
hu et al.
<cit>  and borozan et al.
<cit>  proposed methods that extend traditional effect-size models to combine information from different microarray studies as a way to evaluate or unify lists of genes differentially expressed across them.
another approach combines data from multiple microarray studies  in a melting pot of data and analyzes them as a single data set.
kim and webster  <cit>  used public databases containing microarray data and biological traits on cytoarchitectural abnormalities from the same samples of patients belonging to three groups of major mental disorders plus a control group.
their study used gene expression data measured through two array types, the affymetrix human genome u <dig> set a and the affymetrix human genome 95av <dig>  and the authors carried out a correlation analysis between each gene expression and the biological traits of each subject; although not fully described in the paper, it seems the correlation analysis was performed on the pooled data set from independent studies.
subsequent gene ontology   <cit>  enrichment analysis revealed significant overrepresentation of biological processes, such as cellular metabolism, central nervous system development, cell motility, and programmed cell death, in groups of genes that were significantly correlated with biological traits.
mentzen and wurtele  <cit>  and horan et al.
<cit>  have created co-expression networks for arabidopsis thaliana based on parameters of co-expression similarity that were estimated from a large pool of microarray data downloaded from public repositories.
mentzen and wurtele  <cit>  pooled data from  <dig> affymetrix gene chips, distributed across  <dig> independent studies encompassing diverse organs, conditions, and genotypes “to quest the transcriptome in response to a wide variety of environmentally, genetically, and developmentally induced perturbations”.
horan et al.
<cit>  pooled data corresponding to  <dig> affymetrix microarrays divided among  <dig> independent studies.
both works used cluster analysis based on pearson correlation coefficients as a measure of similarity of gene expression profiles from arabidopsis.
mentzen and wurtele  <cit>  analyzed data from  <dig>  gene probes on the gene chip and identified clusters of co-expressed genes as regulons.
horan et al.
<cit>  used clusters to identify groups of co-regulated protein of unknown function and protein of known function encoding genes from arabidopsis.
in both works, go enrichment analysis showed that networks based on gene-to-gene correlations estimated from pooling data from multiple microarray studies were not random.
a similar approach has been used to obtain regulon information from a human transcriptomic network derived from almost  <dig>  microarrays  <cit> .
this analysis also showed a non random functional distribution of regulons.
from a statistical standpoint, combining data from independent microarray studies into a large pool as a single set can be acceptable if data homogeneity can be ensured across studies.
yet, this condition is nearly impossible to ensure considering that significant data heterogeneity is reported even for completely replicated microarray experiments carried out by the same lab  <cit> .
nevertheless, it can be argued that go enrichment implies meaningful biology and significant go enrichment has been shown for networks created from pooled data  <cit> .
moreover, information gathered through a single data set analysis has led to gene function knowledge discovery  <cit> .
the objective of this study was to perform a comprehensive analysis of pearson correlation coefficients estimated from pooling heterogeneous groups of data  in a large-scale gene expression analysis of publicly available affymetrix microarrays and compare it to the analysis  that combines statistical results of individual groups .
our study included two specific objectives:  to determine the specific effect of different means and heteroskedasticity across the many groups comprising a pool of microarray data on the sign and magnitude of gene-to-gene pearson correlation coefficients obtained from the pool of data, and  to quantify the extent of bias in gene-to-gene pearson correlation coefficients obtained from a pool of heterogeneous groups of microarray data.
in the “methods” section of this article, we describe the statistical theory that we applied to analyze the components of pearson correlation coefficients obtained from a pool of heterogeneous microarray groups.
the “simulation study” section provides results of a study that further tests the specific effect on pearson correlation coefficients of only mean differences, and only heteroskedasticity across n >  <dig> groups, and the validity of our methodology when groups have a small number of elements.
in the section “application to experimental microarray data” we illustrate the results predicted by both theory and simulation with data from  <dig> microarray experiments.
at the end of this section, we provide an assessment of the bias of correlation coefficients estimated across a pool of heterogeneous groups of microarray data.
we discuss our results and summarize our conclusions in the last section.
methods
dissecting components of the pearson correlation coefficient obtained from a pool of microarray data
hassler and thadewald  <cit>  developed the asymptotic formulation to quantify and explain differences between the pearson correlation coefficient estimated from combining two heterogeneous groups into one pool and the pearson correlation coefficients estimated within each group.
they illustrated the problem with a set of measurements on height  and weight  reported by  <dig> first-year college students with a roughly even number of males and females.
as the authors emphasized in their work, male and female groups were not homogeneous because “male students are taller and heavier than female students, and variation around the mean also differs between the groups”.
the generalization of hassler and thadewald’s  <cit>  asymptotic analysis for n heterogeneous groups  is provided in equation  <dig>  for the purpose of applying their theoretical work to analyze correlation coefficients obtained from a pool of microarray data, we consider n heterogeneous groups of gene expression data measured through microarrays.
each group of data can be described as a matrix mi of g genes by ni columns .
we assume that expression levels of any given gene pair xy within each group, i.e.
x,y ∈ mi, are bivariate random normal variables that are identically distributed with means μxy,i =  and variance-covariance matrix Σxy,i=σx,i2σxy,iσxy,iσy,i <dig>  ∀ i =  <dig> n. therefore, heterogeneities across n groups of microarray data are characterized by μxy,i ≠ μxy,j and/or Σxy,i ≠ Σxy,i, for i ≠ j.
the limit in probability of the pearson correlation coefficient between expressions of genes x and y, rxy, obtained from a pool of n heterogeneous groups, as ni → ∞, is given by the expression in equation 1:
rxy→pτxy=∑i=1nλiσxy,i+∑i=1n∑j=i+1nλiλjδxδy
where λi=ni∑i=1nni represents the weight of the number of microarrays of each group and the terms δx <dig> and δy <dig> correspond to the average of expression level variances weighted by ni plus the weighted average of the square of mean differences across n groups for genes x and y, respectively .
hence, the limit in probability of gene-to-gene pearson correlation coefficients obtained from combining heterogeneous groups of microarray data is a mixture of the weighted average of all covariances across n groups plus the weighted average of the cross product of the mean differences of genes x and y across n groups.
both terms are divided by a combination of the average of variances of genes x and y and the mean differences of genes x and y across n groups.
δx2=∑i=1nλiσx,i2+∑i=1n∑j=i+1nλiλj <dig>
δy2=∑i=1nλiσy,i2+∑i=1n∑j=i+1nλiλj <dig>
RESULTS
simulation study
this section presents results of a study using simulated data that had the purpose of further investigating correlation coefficients obtained from a pool of n groups under the following specific conditions:  occurrence of only mean differences across n >  <dig> groups and  occurrence of only heteroskedasticity across n >  <dig> groups.
our study also evaluated the validity of using estimates of the asymptotic expression given in equation  <dig> to explain the components of the pooled correlation coefficient when groups comprising the pool of data have a small number of elements ni.
we performed simulations in r  <cit>  for a generic gene pair xy following the procedure detailed in 1–4:
<dig>  for each simulated group i, i =  <dig> n, generate ni data pairs from a multivariate normal distribution with parameters μxy,i and Σxy,i .
each simulated group is a matrix of  <dig> rows  by ni columns .
<dig>  combine simulated data of n groups into one pool of  <dig> rows  by ∑i=1nni columns.
<dig>  obtain the pearson correlation coefficient from the pool of data .
<dig>  repeat steps 1– <dig> above  <dig> times; results are presented as averages over  <dig> repetitions.
as a control, we first performed an experiment with parameters μxy,i =  and Σxy,i=1001for all groups i =  <dig> .
this simulation provided nearly zero correlation coefficients , thus reassuring us that our simulation procedure worked as expected.
simulation of only mean differences across a pool of n groups
first, we analyzed the case in which heterogeneities across n groups of simulated data were due only to mean differences, i.e.
μxy,i ≠ μxy,j for i ≠ j, but variance-covariances remained constant, i.e.
Σxy,i = Σxy.
we first simulated the case of zero correlation within each group and the effect of differing means  in only two of the n groups.
the simulation results for the set of parameters μxy,1 = , μxy,2 = , −10 ≤ α ≤  <dig>  and μxy,i =  for i ≥ 3; Σxy,i= <dig>  ni =  <dig>  λi = λ,  <dig>  ≤ λ ≤  <dig> , ∀ i =  <dig> n, and 10 ≤ n ≤  <dig> are shown in figure  <dig>  the pooled correlation coefficients rxy shown in figure  <dig> are positive for α <  <dig>  negative for α ≥  <dig>  and nearly zero for α =  <dig>  even though data pairs in each group were drawn from populations with zero correlations, rxy ranged from − <dig>  to  <dig> .
also shown in figure  <dig> is the non-linear relationship between rxy and α as well as between rxy and λ. not surprisingly, coefficients rxy increase as λ increases.
the smooth curves observed in figure  <dig>  are explained by the asymptotic formulation of equation  <dig> as written for the set of population parameters used in this simulation case study :
τxy=−λ2α2−4λλ2α2+λ2+4λ+ <dig>
through equation  <dig> one can see that −4λ24λ2+4λ+1≈ <dig> for α =  <dig> and  <dig>  ≤ λ ≤  <dig> ; τxy >  <dig> for α <  <dig> because the term −4λ >  <dig> and dominates the term −λ2α2 < 0; τxy <  <dig> for α >  <dig> because −4λ <  <dig> and −λ2α2 <  <dig>
secondly, we simulated the case in which the means of a gene pair xy differ for all n groups but the correlation ρxy ≠  <dig> assumes the same value within each group.
the simulation results for the set of parameters μxy,i = ), Σxy,i=1ρxyρxy <dig>  − <dig>  ≤ ρxy ≤  <dig> , ni =  <dig>  λi = λ,  <dig>  ≤ λ ≤  <dig> , for i =  <dig> n, and 10 ≤ n ≤  <dig> are shown in figure 2a.
the correlation coefficients rxy shown in figure 2a are always negative, ranging from − <dig>  < rxy < − <dig> , even when the true correlation within groups ρxy was positive.
the asymptotic formulation of equation  <dig> as written for the set of parameters used in this simulation case study is given in equation 5:
τxy=ρ−λ2∑i=1n∑j=i+1n21+λ2∑i=1n∑j=i+1n <dig>
through equation  <dig> one can see that τxy assumes values of nearly − <dig> for 10 ≤ n ≤  <dig>  a visualization of the problem considered in our second simulation case study is shown  through a scatterplot of expression data of a gene pair xy simulated for  <dig> groups according to the parameters μxy,i = ), ρxy =  <dig> , and ni =  <dig>  for i =  <dig> .
as shown in figure 2b, even though there is a positive trend between the expressions of gene x and gene y within each of the  <dig> groups, the trend between expression of the gene pair xy across the pool of all groups is negative.
simulation of only heteroskedasticity across a pool of n groups
through simulation, we analyzed the effect of variations in Σxy,i=σx,i2σxy,iσxy,iσy,i2among n groups, keeping μxy,i = μ constant.
simulation of the case where σxy,i =  <dig> but σx,i2≠σx,j <dig> and σy,i2≠σy,j <dig> for i ≠ j resulted in nearly zero correlation coefficients from the pool of n groups .
this result agrees with equation  <dig>  which predicts zero correlation if covariances and mean differences of a gene pair across all groups are zero, even when variances differ across groups.
we performed a simulation experiment in which the variance of gene x changes across n groups, but the variance of gene y and the correlation between genes x and y remain constant across n groups.
results of this experiment for the set of parameters μxy,i = , Σxy,i=σx,i2ρσx,iρσx,i <dig>  σx,i2=i <dig>  − <dig>  ≤ ρ ≤  <dig> , ni =  <dig>  λi = λ, and  <dig>  ≤ λ ≤  <dig>  for i =  <dig> n and 10 ≤ n ≤  <dig> are shown in figure  <dig>  the range of the pooled correlations is − <dig>  ≤ rxy ≤  <dig>  , whereas the range of the true correlations within groups is − <dig>  ≤ ρ ≤  <dig>  .
figure  <dig> shows a clear linear relationship between rxy and ρ, in which |rxy| < |ρ|.
the slope between rxy and ρ, estimated through ordinary least squares, is  <dig> .
equation  <dig> gives the asymptotic formulation for the set of parameters used in this simulation case study:
τxy=ρσx¯σx2¯
where σ¯x=∑i=1nσx,in and σx2¯=∑i=1nσx,i2n.
equation  <dig> also shows a linear relationship between τxy and ρ with a slope of σx¯σx2¯= <dig> , for σx,i = i <dig>  i =  <dig> n, and 10 ≤ n ≤  <dig>
there is also a linear relationship between the average of correlation coefficients within each group r¯xy and ρ , in which − <dig> <r¯xy<+ <dig> .
in addition, as shown in figure  <dig>  |r¯xy|>|rxy|.
hence, it can be easily inferred that the mean squared error between rxy and ρ  is greater than the mean squared error between r¯xy and ρ. therefore, pools of n simulated groups marked by only heteroskedasticity provide less efficient estimates of the correlation across groups than combining the n groups’ correlation coefficients into an average.
we analyzed the mean squared error between τxy and τ^xy versus ni, the number of simulated elements in each group, for 10 ≤ ni ≤  <dig> and n =  <dig> .
τxy was obtained from plugging population parameters μxy,i and Σxy,i into equation  <dig>  whereas τ^xy was based on a combination of parameters of each group μ^xy,i and Σ^xy,i .
the mean squared error ranged from  <dig>  for ni =  <dig> to  <dig>  for ni =  <dig>  the correspondence between τxy and τ^xy was good even for ni =  <dig>  a small number of simulated elements per group.
application to experimental microarray data
our simulation study showed that pearson correlation coefficients obtained from a pool of data coming from groups that have different means are explained solely by mean differences across groups.
furthermore, we showed that pooling data marked by only heteroskedasticity provides less efficient estimates of correlation coefficients than does a classical meta-analysis approach of combining correlation coefficients into an average.
the following analysis of experimental microarray data illustrates the results predicted by both theory and simulation.
example data set
the example data set of this work includes the raw expression data from  <dig> affymetrix ath <dig> gene chips  from atgenexpress  <cit> .
cel files are also available from the arabidopsis information resourse   <cit> ; see table  <dig> for the experiment’s id on tair and additional file  <dig> for detailed information about treatment conditions and number of biological replicates.
these data come from  <dig> experiments that explored the effect of  <dig> types of abiotic stress on rna accumulation in shoot and root of  <dig> day-old arabidopsis thaliana seedlings .
experiments followed a 3-factorial design with treatment , plant material , and time post-treatment as factors  <cit> .
seven different research groups located at different institutions across germany performed experiments; microarray data were generated at the german resource center for genome research  .
each microarray experiment is described by its tair  <cit>  id and type of abiotic stress; cel files of each experiment can be located on tair  <cit>  through its id.
total number of downloaded cel files from each experiment is shown on column cel files.
the combination abiotic stress/plant part gives the  <dig> groups used for the asymptotic analysis of pearson correlation coefficients.
the column ni shows number of elements in each of  <dig> groups comprising the large expression matrix.
the columnni∗ shows number of elements in each of  <dig> groups comprising the large matrix of residuals.
experimental data analysis
we imported data from cel files into the r environment  <cit>  and processed the data with mas <dig> from the open-source bioconductor r package affy  <cit> .
following the methodology described in horan et al.
<cit> , we did not screen our example data set for quality of biological replicates, and therefore no outliers were removed.
we followed this procedure because the same data from the  <dig> experiments of atgenexpress  <cit>  were also part of the larger data set used in the work by horan et al.
<cit> .
as described in the methodology of mentzen and wurtele  <cit> , all data were subsequently normalized using the median absolute deviation method as performed by the function normalizebetweenarrays  from the open-source bioconductor r package limma  <cit> .
we obtained mean values of biological replicates after a log transformation  of the normalized expression data.
because the two treatment conditions “genotoxic stress applied to root  <dig> hour post-treatment” and “heat control applied to shoot  <dig> hours post-treatment” had data for only one biological replicate, their expression measurements were used as mean values .
thereafter, mean values of biological replicates were combined into one large expression matrix  encompassing  <dig> columns and  <dig>  rows .
all but two columns of the large expression matrix resulted from averaging data of two or three biological replicates .
gene-to-gene pearson product-moment correlation coefficients  were obtained from the large expression matrix  with the r function cor.
<cit> .
mean differences across a pool of microarray data
we used estimates of the asymptotic expression given in equation  <dig> to examine the makeup of pearson correlation coefficients obtained from pooling the means of biological replicates of different experimental conditions into one large expression matrix.
in order to accomplish this task, we classified data in columns of the large expression matrix into  <dig> groups.
each group had gene expression values from either root or shoot in each of nine types of abiotic stress treatments .
we adopted this procedure because an exploratory analysis showed clear mean differences in gene levels expressed in root or shoot.
the light stress experiment was for entire seedlings, and our analysis did not show mean differences that would justify further division of the data from this experiment.
each group’s name and its corresponding number of elements ni, for groups i =  <dig> , are given in table  <dig> .
data across  <dig> groups were obviously not homogeneous because each group corresponds to a combination of the type of abiotic stress and the plant material, which surely would have an effect on the total group mean of a gene expression.
in addition, data within groups cannot strictly be considered homogeneous either because gene mean expression values within groups correspond to different time points post application of abiotic stress/control treatments .
because our exploratory analysis indicated that the grand mean expression level of genes within groups seemed to dominate over means of all other treatment effects , we considered data within groups as roughly homogeneous.
we used the procedure described in steps  <dig> through  <dig> below to make a diagnostic of rxy obtained from a pool of gene expression data coming from  <dig> heterogeneous groups, where μxy,i ≠ μxy,j and/or Σxy,i ≠ Σxy,j, for i ≠ j.
<dig>  for a given gene pair xy, obtain estimatesμ^xy,i= andΣ^xy,i=sx,i2sxy,isxy,isy,i <dig> for all groups i =  <dig> .
herex¯i andy¯i are, respectively, group means of expression levels of genes x and y, andsx,i <dig> sy,i <dig>  and sxy,i are group variances and covariances of expression levels of genes x and y, respectively.
<dig>  estimate asymptotic coefficientsτ^xy as
τ^xy=s¯xy+dxydxdy
where
s¯xy=∑i=119λisxy,i
dxy=∑i=119∑j>i19λiλj
dx2=sx2¯+∑i=119∑j>i19λiλj <dig>
dy2=sy2¯+∑i=119∑j>i19λiλj <dig>
sx2¯=∑i=119λisx,i <dig>
sy2¯=∑i=119λisy,i <dig>
<dig>  use the residual error  to compare pearson correlation coefficients as obtained from the large expression matrix and coefficients as estimated through equation  <dig>  which are based on parameter estimates of  <dig> groups of data.
<dig>  small residual errors indicate good agreement between rxy andτ^xy.
as a result, the two components of equation  <dig>  can explain the signs and magnitudes of rxy, the pearson correlation coefficient obtained from the large expression matrix.
pearson correlation coefficients obtained from the pooled expression data
we first obtained pearson correlation coefficients for pairwise combinations of all  <dig>  genes present in the large expression matrix, which resulted in more than  <dig> million coefficients.
the pearson correlation coefficients ranged from − <dig>  to  <dig>  with  <dig>  as the mean value, and coefficients showed a symmetric distribution around zero; roughly 10% of these coefficients were greater than  <dig>  or less than − <dig>  .
because all genes present in the large expression matrix provide more than  <dig> million pairwise correlation coefficients, we used a subset of  <dig> randomly selected genes and all their  <dig>  pairwise correlation coefficients to illustrate potential problems with gene pairwise correlation coefficients estimated from a pool of microarray data.
pearson correlation coefficients from the pooled expression data  ranged from − <dig>  to  <dig>  with  <dig>  as the mean value, and coefficients showed a symmetric distribution around zero; roughly 10% of these coefficients were greater than  <dig>  or less than − <dig> , as shown in the histogram of figure  <dig>  the asymptotic coefficientsτ^xy, estimated according to equation  <dig>  ranged from − <dig>  to  <dig>  with  <dig>  as the mean value .
the histogram of residual errors  shows a bimodal distribution in which the mean value of negative residual errors is − <dig>  and the mean value of positive residual errors is  <dig> .
the bimodality of the residual errors implies that|rxy|>|τ^xy|.
in addition, the plot of <dig> versus rxy  shows that residual errors are smaller closer to extreme values and reach a maximum around ± <dig> .
the bimodality observed in figure 5a and the shape observed in figure 5b closely follow the bimodality and shape of the bias between the pearson estimator and the true correlation of a population ρ, which is approximately ρ/  <cit> .
the bias ρ/ is maximized as ρ assumes a value around ± <dig> .
the analysis of all  <dig>  pairwise correlations of  <dig> randomly selected genes revealed good agreement between rxy andτ^xy, despite the approximations we made about homogeneity of data within groups and the relatively low number of elements in each group; in our example data set 12 ≤ ni ≤  <dig>  whereas hassler and thadewald’s example data set had around  <dig> elements in each of two groups  <cit> .
therefore, our analysis reassured us that the pearson correlation coefficients obtained from the large expression matrix can be explained by heterogeneities due to different means and variances-covariances across the  <dig> groups we used to classify our example data set.
next, we show the influence of each term of equation  <dig> on signs and magnitudes of rxy, the pearson correlation coefficients obtained from the large expression matrix.
the plot of rxy versuss¯xy  shows that rxy ranges from − <dig> to + <dig> for negative and positive values ofs¯xy.
therefore, positive or negative covariances of a gene pair within each of the  <dig> groups have no effect on positive or negative correlations estimated from the large expression matrix.
conversely, the “s” shape observed in plot of rxy versus dxy  indicates that positive or negative mean-differences dxy  of a gene pair across the  <dig> groups are the sole determinant of the sign of rxy, i.e.
dxy > 0 ⇒ rxy >  <dig> and dxy < 0 ⇒ rxy <  <dig>  the magnitude of rxy is due mostly to mean differences because the correlation between rxy anddxydxdy is  <dig>  , whereas the correlation between rxy ands¯xydxdy is  <dig>  .
because |rxy| >  <dig>  obtained from pools of microarray data has been used as the cut-off value representing a strong association between gene pairs  <cit> , we computed the percentage contribution of the covariance and mean differences terms on the magnitude of |rxy| ≥  <dig> , i.e.s¯xyrxydxdy+dxyrxydxdy≈ <dig>  there were  <dig>  correlation coefficients with roughly equal numbers distributed in the rxy < − <dig>  and rxy >  <dig>  categories.
the median ofs¯xyrxydxdy% was  <dig> % with 50% of the data showing values between  <dig> % and  <dig> %.
conversely, the median ofdxyrxydxdy% was  <dig> % with 50% of the data showing values between  <dig> % and  <dig> %.
a combination of correlation coefficients between expression profiles within each group, given byr¯xy=∑i=119λirxy,i, ranged from− <dig> <r¯xy< <dig>  with  <dig>  as the mean value.
a direct comparison between rxy andr¯xy showed a correlation coefficient of  <dig> .
by applying the asymptotic theory developed by hassler and thadewald  <cit>  to the pearson correlation coefficients obtained from the large expression matrix, we showed that differences in means across  <dig> heterogeneous groups of data is the main factor determining the magnitude and sign of coefficients of  <dig>  gene pairs.
as previously shown by hassler and thadewald  <cit> , this result corroborates that gene pairwise correlation coefficients estimated from a pool of microarray data do not measure “the closeness of linear relationship”  <cit>   between expressions of a gene pair.
instead, they measure the extent of mean differences of a gene pair across different groups comprising the pool of data.
heteroskedasticity across a pool of microarray data
here, we examine the case in which pearson correlation coefficients are obtained from a pool of microarray data in which only gene pairwise variances-covariances differ across groups of data, i.e.
Σxy,i ≠ Σxy,j for i ≠ j. in this case, equation  <dig> can be written as
rxy∗→pτxy∗=∑i=1nλiσxy,i∑i=1nλiσx,i2∑i=12λiσy,i <dig>
for instance, heteroskedasticity could occur in a situation in which data from completely replicated microarray experiments are pooled to be examined as one data set.
as was reported in the work of goldstein et al.
<cit> , data variability could differ substantially across replicated microarray experiments.
in order to attain only heteroskedasticity across the  <dig> groups of our example data set, we removed the effect of varied experimental conditions on expressions of genes within each group.
for this purpose, we fitted linear models to genes  and obtained their residuals.
following the methodology for differential expression of genes proposed by smith  <cit> , we modeled the expression level of all genes in group i, here represented by a matrix yi, with a systematic treatment effect  plus error, i.e.
yi=ziβi+εi
for i =  <dig> .
we assumed that εi∼n, where Σi is the variance-covariance matrix of all genes in each group i =  <dig> .
we obtained residuals as
ε^i=yi−ziβ^i
whereβ^i was estimated using the open-source bioconductor r package limma  <cit> .
this approach is equivalent to subtracting expression levels of each biological replicate from their mean values.
we used linear models because they are well known by the community who works with differential expression of microarrays measurements and the process of obtaining their residuals is easy and automatic through the use of the limma package  <cit> .
we combined all gene expression residuals from the  <dig> groups into one pool of residuals .
expression levels of the two treatment conditions “genotoxic stress applied to root  <dig> hour post-treatment” and “heat control applied to shoot  <dig> hours post-treatment” could not be used in the analysis of residuals because they had only one biological replicate .
this explains why the matrix of residuals has  <dig> columns instead of  <dig> columns.
we repeated the analysis described in steps 1– <dig> from the section “application to experimental microarray data” for the data in the large matrix of residuals.
pearson correlation coefficients estimated from the pooled residuals
here we show results of the analysis involving the large matrix of residuals  for the same subset of  <dig> genes used in the analysis of the large expression matrix.
pearson correlation coefficients  of all  <dig>   <dig> pairwise combinations of  <dig> genes obtained from the large matrix of residuals ranged from − <dig>  to  <dig>  with  <dig>  as the mean value.
their asymptotic counterparts , estimated according to equation  <dig>  ranged from − <dig>  to  <dig>  with  <dig>  as the mean value.
the combination of covariances within each of the  <dig> groups, i.e.∑i=119λisxy,i, determined the sign ofrxy∗ because all pairwise mean differences among groups were zero .
we then comparedrxy∗≈∑i=119λisxy,i∑i=119λisx,i2∑i=119λisy,i <dig> to the weighted average of correlations obtained within each of the  <dig> groups of residuals, i.e.r¯xy∗=∑i=119λirxy,i∗;r¯xy∗ ranged from − <dig>  to  <dig>  with  <dig>  as the mean value.
in addition, we observed a strong linear relationship betweenrxy∗ estimated from the large matrix of residuals andr¯xy∗, with a correlation equal to  <dig> .
therefore, the pearson correlation coefficients obtained from the large matrix of residuals  also measure a linear relationship between the expression residuals of a gene pair.
bias of correlation coefficients obtained across  <dig> groups of microarray data
we provide here a performance metric for the correlation coefficients estimated across the  <dig> groups of microarray data by assessing their bias with respect to coefficients within each of the  <dig> groups.
we quantified bias as in equation 15:
b=∑i=119λi <dig>
whereρ^xy represents the correlation point estimate of a gene pair xy across the  <dig> groups of microarray data andρ^xy,i represents its counterparts within each group.
we evaluated the bias  of each of the  <dig>  gene pairs’ correlation coefficients that were obtained according to: ρ^xy=rxy, the pearson correlation coefficients obtained directly from the large expression matrix ; ρ^xy=r¯xy, the average of correlations between expression profiles within i =  <dig>  groups comprising the large expression matrix; ρ^xy=rxy∗, the pearson correlation coefficients estimated directly from the large matrix of residuals ; and ρ^xy=r¯xy∗, the average of correlations between expression residuals within i =  <dig>  groups comprising the large matrix of residuals.
for the large expression matrix,ρ^xy,i=rxy,i is the pearson correlation coefficient between expression profiles within each of  <dig> groups comprising the large expression matrix, whereas for the large matrix of residuals,ρ^xy,i=rxy,i∗ is the pearson correlation coefficient between expression residuals within each of  <dig> groups comprising the large matrix of residuals.
table  <dig> gives the statistical summaries of the values obtained for b,b,b, andb.
b – bias of the pearson correlation coefficients estimated directly from the large expression matrix;b – bias of the average of correlations between expression profiles within i =  <dig>  groups; b – bias of the pearson correlation coefficients estimated directly from the pooled residuals;b – bias of the average of correlations between expression residuals within i =  <dig>  groups.
the analysis involving the data in the large expression matrix  resulted in consistently larger statistical summaries of b than did those ofb.
in addition, the maximum value of b is twice as much the maximum value ofb .
for the large matrix of residuals , the values ofb shown in table  <dig> are slightly larger than are the values ofb.
moreover, more information can be grasped through the visualization of biases versus coefficients, as shown in figures 7a and 7b.
the trend shown in the plot of b versus rxy , where b increases as |rxy| approaches ± <dig>  is very distinct from that shown in the plot of b versusr¯xy , where b decreases asr¯xy approaches ± <dig>  indeed, the mean value b for |rxy| >  <dig>  is  <dig> , whereas the mean value ofb for|r¯xy|> <dig>  is  <dig> .
the visualization of biases involving the large matrix of residuals showed a roughly random pattern betweenb and|rxy∗|, as|rxy∗| decreases to zero .
the plot ofb versusr¯xy∗ showed a pattern similar to the one observed in figure 7b, whereb decreases asr¯xy∗ approaches ± <dig> and bothb,b show maximum values around zero .
the plot ofb versusb  reveals all data points lying below the diagonal, thus implying thatb<b, ∀r¯xy∗,rxy∗.
this result corroborates that, in the case of only heteroskedasticity across the  <dig> groups of microarray data, the combination of correlation coefficients performs better than pooling data.
discussion and
CONCLUSIONS
discussion
in this study, we performed a comprehensive analysis of pearson correlation coefficients obtained from combining data of  <dig> heterogeneous groups of experimental microarray data into one large pool.
by applying the theory developed by hassler and thadewald  <cit>  to our example data set, we determined the specific effect of mean differences and heteroskedasticity across the  <dig> groups on the sign and magnitude of the pooled coefficients.
in addition, we provided a performance metric for correlation coefficients by quantifying their biases.
we quantified the bias of the correlation coefficient of a gene pair through the mean squared error between its estimate across a pool of groups and its estimates within groups.
a similar method has been used by hunter and schmidt to assess the variance of their meta-analysis estimator of the pearson correlation coefficient across independent studies  <cit> .
we evaluated the bias of gene-to-gene correlations estimated according to the following two methods:  by combining  <dig> groups of microarray data into a large pool to be analyzed as a single data set  and  by combining correlation coefficients of each of  <dig> groups of microarray data into an average weighted by the number of elements in each group, which corresponds to the hunter-schmidt meta-analysis estimator of the pearson correlation coefficient across independent studies  <cit> .
the data used in this study came from  <dig> microarray experiments  carried out by seven different laboratories distributed across germany that followed the same experimental protocol; these are a subset of the large pool of microarray data found in the study of horan et al.
<cit> .
experiments followed a 3-factorial design with treatment , tissue , and time post-treatment as factors  <cit> .
mean differences within and across experiments were a matter-of-fact because statistically significant differences in gene expression of several types of abiotic stress versus control treatment were reported in kilian et al.’s study  <cit> .
we expected differences due to variability across experiments to remain after removing mean differences because of reported difficulties in the reproduction of microarray studies  <cit> .
therefore, homogeneity cannot be ensured across experiments, and combining the means or residuals of biological replicates of the  <dig> experiments into a large pool as a single set is not sound from a statistical viewpoint.
the analysis of the components of the correlation coefficients obtained from the large expression matrix corroborated the results predicted by both theory and simulation that variances-covariances within the  <dig> groups had negligible impact on correlation coefficients, but different means across the  <dig> groups had a decisive effect on the sign as well as on the magnitude of coefficients.
coefficients that were greater than  <dig>  or less than − <dig>  showed the largest range of bias .
therefore, large values of the pooled coefficients were an artifact in the sense that they did not communicate a real linear association between the expression profiles of two genes; rather, they appeared because the data were combined into a large pool.
for this reason, large values of the pooled coefficients are in fact an ecological fallacy  <cit> .
we also showed through monte carlo simulation that the structure of different means across a pool of 10 ≤ n ≤  <dig> groups could generate simpson’s paradox.
in our case study simulation shown in figure  <dig>  we showed that even though the correlation within each group was + <dig> , a pool of n  groups provided negative correlation coefficients because the combination of all pairwise mean differences had a negative sign and greater magnitude than the positive covariance of the data within groups.
hassler and thadewald  <cit>  studied simpson’s paradox through the analytical solution of equation  <dig> for n =  <dig>  and showed that the occurrence of mean differences with opposite signs in both correlated variables is a condition for contradictory results between a correlation coefficient that is estimated across or within each of two groups.
we combined residuals from fitting linear models of every gene into a large matrix of residuals .
here we departed from the assumption of independence  and considered a multivariate normal distribution for residuals within groups, with a mean of zero and variance-covariance Σxy,i, i =  <dig> .
the large matrix of residuals gave us the opportunity to evaluate gene pair correlations estimated from a pool of data marked by only heteroskedasticity.
our results showed that correlation coefficients estimated across the  <dig> groups of residuals were closely related to the variance-covariances within groups.
we also found a strong linear relationship between the pearson correlation coefficients obtained from the large matrix of residuals and the coefficients resulting from averaging correlation estimates within groups.
however, the heteroskedasticity of the data in the large matrix of residuals resulted in less efficient estimations of the correlation between a gene pair than did the classical meta-analysis approach of combining correlation coefficients into an average.
these results were corroborated by monte carlo simulations of only heteroskedasticity across n >  <dig> groups of data.
the results shown in this study indicate that the combination of statistical results is best suited for estimating correlations of a gene pair across several microarray studies.
nevertheless, further studies are necessary to assess various methods of combining within-groups gene-to-gene correlation coefficients.

 there is a wealth of information enclosed in the massive amount of microarray data so far accumulated in public repositories. the variety of data sets generated from the assortment of experiments is a major obstacle in the path leading from these data to information. specific issues relating to data heterogeneity across microarray studies include differences across platforms, differences within and between labs, and/or differences among experimental factors such as treatments and tissues. furthermore, concerns regarding integration of studies from multiple sources in general, such as variations in design, research goals, or quality of implementation, add to these issues. inappropriate integration of microarray data from public repositories could lead to equivocal results. the “simpson’s paradox”, which refers to contradictory statistical results obtained when analysis is performed within versus across groups of data, is an example of mishandling of data. blyth gives an example involving the analysis of 2x contingency tables across two groups, and hassler and thadewald also illustrate simpson’s paradox when correlation coefficients are estimated from a pool of two groups versus within each group. in both cases, the paradox can be explained as results are further investigated in light of the specific statistical properties of each group of data. the “ecological fallacy” happens when the correlation of aggregated variables results in a significant relationship that is due only to aggregation rather than to any real association. an early example of an ecological fallacy can be found in gehlke and biehl, whose study of grouping effects in census tract data showed that the magnitude of correlation coefficients of two variables tend to increase as the level of census tract aggregation increases. this problem was later referred to as the “modifiable areal unit problem” and further studied by openshaw and taylor. combining statistical results of independent studies that address similar questions has been a standard procedure in classic meta-analyses. this approach entails analyzing each data set independently and then combining the results, as in a mosaic. meta-analysis of microarray data has been applied in a broader context, as some works include data spanning a wide range of purposes and designs. parmigiani et al., in a quest for a common gene signature across multiple cancer types, developed a statistical method to identify and assess the intersection of multiple gene expression signatures across published cancer-related microarray studies. on the other hand, wirapati et al. and rhodes et al. developed specific meta-analysis methods to integrate gene expression signatures of breast and lung cancer, respectively, across independent studies of microarray data. hu et al. and borozan et al. proposed methods that extend traditional effect-size models to combine information from different microarray studies as a way to evaluate or unify lists of genes differentially expressed across them. another approach combines data from multiple microarray studies in a melting pot of data and analyzes them as a single data set. kim and webster used public databases containing microarray data and biological traits on cytoarchitectural abnormalities from the same samples of patients belonging to three groups of major mental disorders plus a control group. their study used gene expression data measured through two array types, the affymetrix human genome u set a and the affymetrix human genome 95av and the authors carried out a correlation analysis between each gene expression and the biological traits of each subject; although not fully described in the paper, it seems the correlation analysis was performed on the pooled data set from independent studies. subsequent gene ontology enrichment analysis revealed significant overrepresentation of biological processes, such as cellular metabolism, central nervous system development, cell motility, and programmed cell death, in groups of genes that were significantly correlated with biological traits. mentzen and wurtele and horan et al. have created co-expression networks for arabidopsis thaliana based on parameters of co-expression similarity that were estimated from a large pool of microarray data downloaded from public repositories. mentzen and wurtele pooled data from affymetrix gene chips, distributed across independent studies encompassing diverse organs, conditions, and genotypes “to quest the transcriptome in response to a wide variety of environmentally, genetically, and developmentally induced perturbations”. horan et al. pooled data corresponding to affymetrix microarrays divided among independent studies. both works used cluster analysis based on pearson correlation coefficients as a measure of similarity of gene expression profiles from arabidopsis. mentzen and wurtele analyzed data from gene probes on the gene chip and identified clusters of co-expressed genes as regulons. horan et al. used clusters to identify groups of co-regulated protein of unknown function and protein of known function encoding genes from arabidopsis. in both works, go enrichment analysis showed that networks based on gene-to-gene correlations estimated from pooling data from multiple microarray studies were not random. a similar approach has been used to obtain regulon information from a human transcriptomic network derived from almost microarrays. this analysis also showed a non random functional distribution of regulons. from a statistical standpoint, combining data from independent microarray studies into a large pool as a single set can be acceptable if data homogeneity can be ensured across studies. yet, this condition is nearly impossible to ensure considering that significant data heterogeneity is reported even for completely replicated microarray experiments carried out by the same lab. nevertheless, it can be argued that go enrichment implies meaningful biology and significant go enrichment has been shown for networks created from pooled data. moreover, information gathered through a single data set analysis has led to gene function knowledge discovery. the objective of this study was to perform a comprehensive analysis of pearson correlation coefficients estimated from pooling heterogeneous groups of data in a large-scale gene expression analysis of publicly available affymetrix microarrays and compare it to the analysis that combines statistical results of individual groups. our study included two specific objectives: to determine the specific effect of different means and heteroskedasticity across the many groups comprising a pool of microarray data on the sign and magnitude of gene-to-gene pearson correlation coefficients obtained from the pool of data, and to quantify the extent of bias in gene-to-gene pearson correlation coefficients obtained from a pool of heterogeneous groups of microarray data. in the “methods” section of this article, we describe the statistical theory that we applied to analyze the components of pearson correlation coefficients obtained from a pool of heterogeneous microarray groups. the “simulation study” section provides results of a study that further tests the specific effect on pearson correlation coefficients of only mean differences, and only heteroskedasticity across n >  groups, and the validity of our methodology when groups have a small number of elements. in the section “application to experimental microarray data” we illustrate the results predicted by both theory and simulation with data from microarray experiments. at the end of this section, we provide an assessment of the bias of correlation coefficients estimated across a pool of heterogeneous groups of microarray data. we discuss our results and summarize our conclusions in the last section. methods dissecting components of the pearson correlation coefficient obtained from a pool of microarray data hassler and thadewald developed the asymptotic formulation to quantify and explain differences between the pearson correlation coefficient estimated from combining two heterogeneous groups into one pool and the pearson correlation coefficients estimated within each group. they illustrated the problem with a set of measurements on height and weight reported by first-year college students with a roughly even number of males and females. as the authors emphasized in their work, male and female groups were not homogeneous because “male students are taller and heavier than female students, and variation around the mean also differs between the groups”. the generalization of hassler and thadewald’s asymptotic analysis for n heterogeneous groups is provided in equation for the purpose of applying their theoretical work to analyze correlation coefficients obtained from a pool of microarray data, we consider n heterogeneous groups of gene expression data measured through microarrays. each group of data can be described as a matrix mi of g genes by ni columns. we assume that expression levels of any given gene pair xy within each group, i.e. x,y ∈ mi, are bivariate random normal variables that are identically distributed with means μxy,i =  and variance-covariance matrix Σxy,i=σx,i2σxy,iσxy,iσy,i ∀ i =  n. therefore, heterogeneities across n groups of microarray data are characterized by μxy,i ≠ μxy,j and/or Σxy,i ≠ Σxy,i, for i ≠ j. the limit in probability of the pearson correlation coefficient between expressions of genes x and y, rxy, obtained from a pool of n heterogeneous groups, as ni → ∞, is given by the expression in equation 1: rxy→pτxy=∑i=1nλiσxy,i+∑i=1n∑j=i+1nλiλjδxδy where λi=ni∑i=1nni represents the weight of the number of microarrays of each group and the terms δx and δy correspond to the average of expression level variances weighted by ni plus the weighted average of the square of mean differences across n groups for genes x and y, respectively. hence, the limit in probability of gene-to-gene pearson correlation coefficients obtained from combining heterogeneous groups of microarray data is a mixture of the weighted average of all covariances across n groups plus the weighted average of the cross product of the mean differences of genes x and y across n groups. both terms are divided by a combination of the average of variances of genes x and y and the mean differences of genes x and y across n groups. δx2=∑i=1nλiσx,i2+∑i=1n∑j=i+1nλiλj δy2=∑i=1nλiσy,i2+∑i=1n∑j=i+1nλiλj  simulation study this section presents results of a study using simulated data that had the purpose of further investigating correlation coefficients obtained from a pool of n groups under the following specific conditions: occurrence of only mean differences across n >  groups and occurrence of only heteroskedasticity across n >  groups. our study also evaluated the validity of using estimates of the asymptotic expression given in equation to explain the components of the pooled correlation coefficient when groups comprising the pool of data have a small number of elements ni. we performed simulations in r for a generic gene pair xy following the procedure detailed in 1–4: for each simulated group i, i =  n, generate ni data pairs from a multivariate normal distribution with parameters μxy,i and Σxy,i. each simulated group is a matrix of rows by ni columns. combine simulated data of n groups into one pool of rows by ∑i=1nni columns. obtain the pearson correlation coefficient from the pool of data. repeat steps 1– above times; results are presented as averages over repetitions. as a control, we first performed an experiment with parameters μxy,i =  and Σxy,i=1001for all groups i = . this simulation provided nearly zero correlation coefficients, thus reassuring us that our simulation procedure worked as expected. simulation of only mean differences across a pool of n groups first, we analyzed the case in which heterogeneities across n groups of simulated data were due only to mean differences, i.e. μxy,i ≠ μxy,j for i ≠ j, but variance-covariances remained constant, i.e. Σxy,i = Σxy. we first simulated the case of zero correlation within each group and the effect of differing means in only two of the n groups. the simulation results for the set of parameters μxy,1 =, μxy,2 =, −10 ≤ α ≤  and μxy,i =  for i ≥ 3; Σxy,i= ni =  λi = λ,  ≤ λ ≤ , ∀ i =  n, and 10 ≤ n ≤  are shown in figure the pooled correlation coefficients rxy shown in figure are positive for α <  negative for α ≥  and nearly zero for α =  even though data pairs in each group were drawn from populations with zero correlations, rxy ranged from − to. also shown in figure is the non-linear relationship between rxy and α as well as between rxy and λ. not surprisingly, coefficients rxy increase as λ increases. the smooth curves observed in figure are explained by the asymptotic formulation of equation as written for the set of population parameters used in this simulation case study: τxy=−λ2α2−4λλ2α2+λ2+4λ+ through equation one can see that −4λ24λ2+4λ+1≈ for α =  and  ≤ λ ≤ ; τxy >  for α <  because the term −4λ >  and dominates the term −λ2α2 < 0; τxy <  for α >  because −4λ <  and −λ2α2 <  secondly, we simulated the case in which the means of a gene pair xy differ for all n groups but the correlation ρxy ≠  assumes the same value within each group. the simulation results for the set of parameters μxy,i = ), Σxy,i=1ρxyρxy −  ≤ ρxy ≤ , ni =  λi = λ,  ≤ λ ≤ , for i =  n, and 10 ≤ n ≤  are shown in figure 2a. the correlation coefficients rxy shown in figure 2a are always negative, ranging from −  < rxy < −, even when the true correlation within groups ρxy was positive. the asymptotic formulation of equation as written for the set of parameters used in this simulation case study is given in equation 5: τxy=ρ−λ2∑i=1n∑j=i+1n21+λ2∑i=1n∑j=i+1n through equation one can see that τxy assumes values of nearly − for 10 ≤ n ≤  a visualization of the problem considered in our second simulation case study is shown through a scatterplot of expression data of a gene pair xy simulated for groups according to the parameters μxy,i = ), ρxy = , and ni =  for i = . as shown in figure 2b, even though there is a positive trend between the expressions of gene x and gene y within each of the groups, the trend between expression of the gene pair xy across the pool of all groups is negative. simulation of only heteroskedasticity across a pool of n groups through simulation, we analyzed the effect of variations in Σxy,i=σx,i2σxy,iσxy,iσy,i2among n groups, keeping μxy,i = μ constant. simulation of the case where σxy,i =  but σx,i2≠σx,j and σy,i2≠σy,j for i ≠ j resulted in nearly zero correlation coefficients from the pool of n groups. this result agrees with equation which predicts zero correlation if covariances and mean differences of a gene pair across all groups are zero, even when variances differ across groups. we performed a simulation experiment in which the variance of gene x changes across n groups, but the variance of gene y and the correlation between genes x and y remain constant across n groups. results of this experiment for the set of parameters μxy,i =, Σxy,i=σx,i2ρσx,iρσx,i σx,i2=i −  ≤ ρ ≤ , ni =  λi = λ, and  ≤ λ ≤  for i =  n and 10 ≤ n ≤  are shown in figure the range of the pooled correlations is −  ≤ rxy ≤ , whereas the range of the true correlations within groups is −  ≤ ρ ≤ . figure shows a clear linear relationship between rxy and ρ, in which |rxy| < |ρ|. the slope between rxy and ρ, estimated through ordinary least squares, is. equation gives the asymptotic formulation for the set of parameters used in this simulation case study: τxy=ρσx¯σx2¯ where σ¯x=∑i=1nσx,in and σx2¯=∑i=1nσx,i2n. equation also shows a linear relationship between τxy and ρ with a slope of σx¯σx2¯=, for σx,i = i i =  n, and 10 ≤ n ≤  there is also a linear relationship between the average of correlation coefficients within each group r¯xy and ρ, in which − <r¯xy<+. in addition, as shown in figure |r¯xy|>|rxy|. hence, it can be easily inferred that the mean squared error between rxy and ρ is greater than the mean squared error between r¯xy and ρ. therefore, pools of n simulated groups marked by only heteroskedasticity provide less efficient estimates of the correlation across groups than combining the n groups’ correlation coefficients into an average. we analyzed the mean squared error between τxy and τ^xy versus ni, the number of simulated elements in each group, for 10 ≤ ni ≤  and n = . τxy was obtained from plugging population parameters μxy,i and Σxy,i into equation whereas τ^xy was based on a combination of parameters of each group μ^xy,i and Σ^xy,i. the mean squared error ranged from for ni =  to for ni =  the correspondence between τxy and τ^xy was good even for ni =  a small number of simulated elements per group. application to experimental microarray data our simulation study showed that pearson correlation coefficients obtained from a pool of data coming from groups that have different means are explained solely by mean differences across groups. furthermore, we showed that pooling data marked by only heteroskedasticity provides less efficient estimates of correlation coefficients than does a classical meta-analysis approach of combining correlation coefficients into an average. the following analysis of experimental microarray data illustrates the results predicted by both theory and simulation. example data set the example data set of this work includes the raw expression data from affymetrix ath gene chips from atgenexpress. cel files are also available from the arabidopsis information resourse; see table for the experiment’s id on tair and additional file for detailed information about treatment conditions and number of biological replicates. these data come from experiments that explored the effect of types of abiotic stress on rna accumulation in shoot and root of day-old arabidopsis thaliana seedlings. experiments followed a 3-factorial design with treatment, plant material, and time post-treatment as factors. seven different research groups located at different institutions across germany performed experiments; microarray data were generated at the german resource center for genome research. each microarray experiment is described by its tair id and type of abiotic stress; cel files of each experiment can be located on tair through its id. total number of downloaded cel files from each experiment is shown on column cel files. the combination abiotic stress/plant part gives the groups used for the asymptotic analysis of pearson correlation coefficients. the column ni shows number of elements in each of groups comprising the large expression matrix. the columnni∗ shows number of elements in each of groups comprising the large matrix of residuals. experimental data analysis we imported data from cel files into the r environment and processed the data with mas from the open-source bioconductor r package affy. following the methodology described in horan et al., we did not screen our example data set for quality of biological replicates, and therefore no outliers were removed. we followed this procedure because the same data from the experiments of atgenexpress were also part of the larger data set used in the work by horan et al.. as described in the methodology of mentzen and wurtele, all data were subsequently normalized using the median absolute deviation method as performed by the function normalizebetweenarrays from the open-source bioconductor r package limma. we obtained mean values of biological replicates after a log transformation of the normalized expression data. because the two treatment conditions “genotoxic stress applied to root hour post-treatment” and “heat control applied to shoot hours post-treatment” had data for only one biological replicate, their expression measurements were used as mean values. thereafter, mean values of biological replicates were combined into one large expression matrix encompassing columns and rows. all but two columns of the large expression matrix resulted from averaging data of two or three biological replicates. gene-to-gene pearson product-moment correlation coefficients were obtained from the large expression matrix with the r function cor.. mean differences across a pool of microarray data we used estimates of the asymptotic expression given in equation to examine the makeup of pearson correlation coefficients obtained from pooling the means of biological replicates of different experimental conditions into one large expression matrix. in order to accomplish this task, we classified data in columns of the large expression matrix into groups. each group had gene expression values from either root or shoot in each of nine types of abiotic stress treatments. we adopted this procedure because an exploratory analysis showed clear mean differences in gene levels expressed in root or shoot. the light stress experiment was for entire seedlings, and our analysis did not show mean differences that would justify further division of the data from this experiment. each group’s name and its corresponding number of elements ni, for groups i = , are given in table. data across groups were obviously not homogeneous because each group corresponds to a combination of the type of abiotic stress and the plant material, which surely would have an effect on the total group mean of a gene expression. in addition, data within groups cannot strictly be considered homogeneous either because gene mean expression values within groups correspond to different time points post application of abiotic stress/control treatments. because our exploratory analysis indicated that the grand mean expression level of genes within groups seemed to dominate over means of all other treatment effects, we considered data within groups as roughly homogeneous. we used the procedure described in steps through below to make a diagnostic of rxy obtained from a pool of gene expression data coming from heterogeneous groups, where μxy,i ≠ μxy,j and/or Σxy,i ≠ Σxy,j, for i ≠ j. for a given gene pair xy, obtain estimatesμ^xy,i= andΣ^xy,i=sx,i2sxy,isxy,isy,i for all groups i = . herex¯i andy¯i are, respectively, group means of expression levels of genes x and y, andsx,i sy,i and sxy,i are group variances and covariances of expression levels of genes x and y, respectively. estimate asymptotic coefficientsτ^xy as τ^xy=s¯xy+dxydxdy where s¯xy=∑i=119λisxy,i dxy=∑i=119∑j>i19λiλj dx2=sx2¯+∑i=119∑j>i19λiλj dy2=sy2¯+∑i=119∑j>i19λiλj sx2¯=∑i=119λisx,i sy2¯=∑i=119λisy,i use the residual error to compare pearson correlation coefficients as obtained from the large expression matrix and coefficients as estimated through equation which are based on parameter estimates of groups of data. small residual errors indicate good agreement between rxy andτ^xy. as a result, the two components of equation can explain the signs and magnitudes of rxy, the pearson correlation coefficient obtained from the large expression matrix. pearson correlation coefficients obtained from the pooled expression data we first obtained pearson correlation coefficients for pairwise combinations of all genes present in the large expression matrix, which resulted in more than million coefficients. the pearson correlation coefficients ranged from − to with as the mean value, and coefficients showed a symmetric distribution around zero; roughly 10% of these coefficients were greater than or less than −. because all genes present in the large expression matrix provide more than million pairwise correlation coefficients, we used a subset of randomly selected genes and all their pairwise correlation coefficients to illustrate potential problems with gene pairwise correlation coefficients estimated from a pool of microarray data. pearson correlation coefficients from the pooled expression data ranged from − to with as the mean value, and coefficients showed a symmetric distribution around zero; roughly 10% of these coefficients were greater than or less than −, as shown in the histogram of figure the asymptotic coefficientsτ^xy, estimated according to equation ranged from − to with as the mean value. the histogram of residual errors shows a bimodal distribution in which the mean value of negative residual errors is − and the mean value of positive residual errors is. the bimodality of the residual errors implies that|rxy|>|τ^xy|. in addition, the plot of versus rxy shows that residual errors are smaller closer to extreme values and reach a maximum around ±. the bimodality observed in figure 5a and the shape observed in figure 5b closely follow the bimodality and shape of the bias between the pearson estimator and the true correlation of a population ρ, which is approximately ρ/. the bias ρ/ is maximized as ρ assumes a value around ±. the analysis of all pairwise correlations of randomly selected genes revealed good agreement between rxy andτ^xy, despite the approximations we made about homogeneity of data within groups and the relatively low number of elements in each group; in our example data set 12 ≤ ni ≤  whereas hassler and thadewald’s example data set had around elements in each of two groups. therefore, our analysis reassured us that the pearson correlation coefficients obtained from the large expression matrix can be explained by heterogeneities due to different means and variances-covariances across the groups we used to classify our example data set. next, we show the influence of each term of equation on signs and magnitudes of rxy, the pearson correlation coefficients obtained from the large expression matrix. the plot of rxy versuss¯xy shows that rxy ranges from − to + for negative and positive values ofs¯xy. therefore, positive or negative covariances of a gene pair within each of the groups have no effect on positive or negative correlations estimated from the large expression matrix. conversely, the “s” shape observed in plot of rxy versus dxy indicates that positive or negative mean-differences dxy of a gene pair across the groups are the sole determinant of the sign of rxy, i.e. dxy > 0 ⇒ rxy >  and dxy < 0 ⇒ rxy <  the magnitude of rxy is due mostly to mean differences because the correlation between rxy anddxydxdy is, whereas the correlation between rxy ands¯xydxdy is. because |rxy| >  obtained from pools of microarray data has been used as the cut-off value representing a strong association between gene pairs, we computed the percentage contribution of the covariance and mean differences terms on the magnitude of |rxy| ≥ , i.e.s¯xyrxydxdy+dxyrxydxdy≈ there were correlation coefficients with roughly equal numbers distributed in the rxy < − and rxy >  categories. the median ofs¯xyrxydxdy% was % with 50% of the data showing values between % and %. conversely, the median ofdxyrxydxdy% was % with 50% of the data showing values between % and %. a combination of correlation coefficients between expression profiles within each group, given byr¯xy=∑i=119λirxy,i, ranged from− <r¯xy< with as the mean value. a direct comparison between rxy andr¯xy showed a correlation coefficient of. by applying the asymptotic theory developed by hassler and thadewald to the pearson correlation coefficients obtained from the large expression matrix, we showed that differences in means across heterogeneous groups of data is the main factor determining the magnitude and sign of coefficients of gene pairs. as previously shown by hassler and thadewald, this result corroborates that gene pairwise correlation coefficients estimated from a pool of microarray data do not measure “the closeness of linear relationship” between expressions of a gene pair. instead, they measure the extent of mean differences of a gene pair across different groups comprising the pool of data. heteroskedasticity across a pool of microarray data here, we examine the case in which pearson correlation coefficients are obtained from a pool of microarray data in which only gene pairwise variances-covariances differ across groups of data, i.e. Σxy,i ≠ Σxy,j for i ≠ j. in this case, equation can be written as rxy∗→pτxy∗=∑i=1nλiσxy,i∑i=1nλiσx,i2∑i=12λiσy,i for instance, heteroskedasticity could occur in a situation in which data from completely replicated microarray experiments are pooled to be examined as one data set. as was reported in the work of goldstein et al., data variability could differ substantially across replicated microarray experiments. in order to attain only heteroskedasticity across the groups of our example data set, we removed the effect of varied experimental conditions on expressions of genes within each group. for this purpose, we fitted linear models to genes and obtained their residuals. following the methodology for differential expression of genes proposed by smith, we modeled the expression level of all genes in group i, here represented by a matrix yi, with a systematic treatment effect plus error, i.e. yi=ziβi+εi for i = . we assumed that εi∼n, where Σi is the variance-covariance matrix of all genes in each group i = . we obtained residuals as ε^i=yi−ziβ^i whereβ^i was estimated using the open-source bioconductor r package limma. this approach is equivalent to subtracting expression levels of each biological replicate from their mean values. we used linear models because they are well known by the community who works with differential expression of microarrays measurements and the process of obtaining their residuals is easy and automatic through the use of the limma package. we combined all gene expression residuals from the groups into one pool of residuals. expression levels of the two treatment conditions “genotoxic stress applied to root hour post-treatment” and “heat control applied to shoot hours post-treatment” could not be used in the analysis of residuals because they had only one biological replicate. this explains why the matrix of residuals has columns instead of columns. we repeated the analysis described in steps 1– from the section “application to experimental microarray data” for the data in the large matrix of residuals. pearson correlation coefficients estimated from the pooled residuals here we show results of the analysis involving the large matrix of residuals for the same subset of genes used in the analysis of the large expression matrix. pearson correlation coefficients of all pairwise combinations of genes obtained from the large matrix of residuals ranged from − to with as the mean value. their asymptotic counterparts, estimated according to equation ranged from − to with as the mean value. the combination of covariances within each of the groups, i.e.∑i=119λisxy,i, determined the sign ofrxy∗ because all pairwise mean differences among groups were zero. we then comparedrxy∗≈∑i=119λisxy,i∑i=119λisx,i2∑i=119λisy,i to the weighted average of correlations obtained within each of the groups of residuals, i.e.r¯xy∗=∑i=119λirxy,i∗;r¯xy∗ ranged from − to with as the mean value. in addition, we observed a strong linear relationship betweenrxy∗ estimated from the large matrix of residuals andr¯xy∗, with a correlation equal to. therefore, the pearson correlation coefficients obtained from the large matrix of residuals also measure a linear relationship between the expression residuals of a gene pair. bias of correlation coefficients obtained across groups of microarray data we provide here a performance metric for the correlation coefficients estimated across the groups of microarray data by assessing their bias with respect to coefficients within each of the groups. we quantified bias as in equation 15: b=∑i=119λi whereρ^xy represents the correlation point estimate of a gene pair xy across the groups of microarray data andρ^xy,i represents its counterparts within each group. we evaluated the bias of each of the gene pairs’ correlation coefficients that were obtained according to: ρ^xy=rxy, the pearson correlation coefficients obtained directly from the large expression matrix; ρ^xy=r¯xy, the average of correlations between expression profiles within i =  groups comprising the large expression matrix; ρ^xy=rxy∗, the pearson correlation coefficients estimated directly from the large matrix of residuals; and ρ^xy=r¯xy∗, the average of correlations between expression residuals within i =  groups comprising the large matrix of residuals. for the large expression matrix,ρ^xy,i=rxy,i is the pearson correlation coefficient between expression profiles within each of groups comprising the large expression matrix, whereas for the large matrix of residuals,ρ^xy,i=rxy,i∗ is the pearson correlation coefficient between expression residuals within each of groups comprising the large matrix of residuals. table gives the statistical summaries of the values obtained for b,b,b, andb. b – bias of the pearson correlation coefficients estimated directly from the large expression matrix;b – bias of the average of correlations between expression profiles within i =  groups; b – bias of the pearson correlation coefficients estimated directly from the pooled residuals;b – bias of the average of correlations between expression residuals within i =  groups. the analysis involving the data in the large expression matrix resulted in consistently larger statistical summaries of b than did those ofb. in addition, the maximum value of b is twice as much the maximum value ofb. for the large matrix of residuals, the values ofb shown in table are slightly larger than are the values ofb. moreover, more information can be grasped through the visualization of biases versus coefficients, as shown in figures 7a and 7b. the trend shown in the plot of b versus rxy, where b increases as |rxy| approaches ± is very distinct from that shown in the plot of b versusr¯xy, where b decreases asr¯xy approaches ± indeed, the mean value b for |rxy| >  is, whereas the mean value ofb for|r¯xy|> is. the visualization of biases involving the large matrix of residuals showed a roughly random pattern betweenb and|rxy∗|, as|rxy∗| decreases to zero. the plot ofb versusr¯xy∗ showed a pattern similar to the one observed in figure 7b, whereb decreases asr¯xy∗ approaches ± and bothb,b show maximum values around zero. the plot ofb versusb reveals all data points lying below the diagonal, thus implying thatb<b, ∀r¯xy∗,rxy∗. this result corroborates that, in the case of only heteroskedasticity across the groups of microarray data, the combination of correlation coefficients performs better than pooling data. discussion and  discussion in this study, we performed a comprehensive analysis of pearson correlation coefficients obtained from combining data of heterogeneous groups of experimental microarray data into one large pool. by applying the theory developed by hassler and thadewald to our example data set, we determined the specific effect of mean differences and heteroskedasticity across the groups on the sign and magnitude of the pooled coefficients. in addition, we provided a performance metric for correlation coefficients by quantifying their biases. we quantified the bias of the correlation coefficient of a gene pair through the mean squared error between its estimate across a pool of groups and its estimates within groups. a similar method has been used by hunter and schmidt to assess the variance of their meta-analysis estimator of the pearson correlation coefficient across independent studies. we evaluated the bias of gene-to-gene correlations estimated according to the following two methods: by combining groups of microarray data into a large pool to be analyzed as a single data set and by combining correlation coefficients of each of groups of microarray data into an average weighted by the number of elements in each group, which corresponds to the hunter-schmidt meta-analysis estimator of the pearson correlation coefficient across independent studies. the data used in this study came from microarray experiments carried out by seven different laboratories distributed across germany that followed the same experimental protocol; these are a subset of the large pool of microarray data found in the study of horan et al.. experiments followed a 3-factorial design with treatment, tissue, and time post-treatment as factors. mean differences within and across experiments were a matter-of-fact because statistically significant differences in gene expression of several types of abiotic stress versus control treatment were reported in kilian et al.’s study. we expected differences due to variability across experiments to remain after removing mean differences because of reported difficulties in the reproduction of microarray studies. therefore, homogeneity cannot be ensured across experiments, and combining the means or residuals of biological replicates of the experiments into a large pool as a single set is not sound from a statistical viewpoint. the analysis of the components of the correlation coefficients obtained from the large expression matrix corroborated the results predicted by both theory and simulation that variances-covariances within the groups had negligible impact on correlation coefficients, but different means across the groups had a decisive effect on the sign as well as on the magnitude of coefficients. coefficients that were greater than or less than − showed the largest range of bias. therefore, large values of the pooled coefficients were an artifact in the sense that they did not communicate a real linear association between the expression profiles of two genes; rather, they appeared because the data were combined into a large pool. for this reason, large values of the pooled coefficients are in fact an ecological fallacy. we also showed through monte carlo simulation that the structure of different means across a pool of 10 ≤ n ≤  groups could generate simpson’s paradox. in our case study simulation shown in figure we showed that even though the correlation within each group was +, a pool of n groups provided negative correlation coefficients because the combination of all pairwise mean differences had a negative sign and greater magnitude than the positive covariance of the data within groups. hassler and thadewald studied simpson’s paradox through the analytical solution of equation for n =  and showed that the occurrence of mean differences with opposite signs in both correlated variables is a condition for contradictory results between a correlation coefficient that is estimated across or within each of two groups. we combined residuals from fitting linear models of every gene into a large matrix of residuals. here we departed from the assumption of independence and considered a multivariate normal distribution for residuals within groups, with a mean of zero and variance-covariance Σxy,i, i = . the large matrix of residuals gave us the opportunity to evaluate gene pair correlations estimated from a pool of data marked by only heteroskedasticity. our results showed that correlation coefficients estimated across the groups of residuals were closely related to the variance-covariances within groups. we also found a strong linear relationship between the pearson correlation coefficients obtained from the large matrix of residuals and the coefficients resulting from averaging correlation estimates within groups. however, the heteroskedasticity of the data in the large matrix of residuals resulted in less efficient estimations of the correlation between a gene pair than did the classical meta-analysis approach of combining correlation coefficients into an average. these results were corroborated by monte carlo simulations of only heteroskedasticity across n >  groups of data. the results shown in this study indicate that the combination of statistical results is best suited for estimating correlations of a gene pair across several microarray studies. nevertheless, further studies are necessary to assess various methods of combining within-groups gene-to-gene correlation coefficients. 
mean differences across the groups were the main factor determining the magnitude and sign of the pooled coefficients, which showed largest values of bias as they approached ± only heteroskedasticity across the pool of groups resulted in less efficient estimations of correlations than did a classical meta-analysis approach of combining correlation coefficients. the synthesis of information across microarray studies has been performed by combining statistical results of individual studies, or by combining data from multiple studies into a large pool to be analyzed as a single data set. specific issues relating to data heterogeneity across microarray studies, such as differences within and between labs or differences among experimental conditions, could lead to equivocal results in a melting pot approach. the combination of statistical results is best suited for synthesizing the correlation between expression profiles of a gene pair across several microarray studies. these results were corroborated by simulation studies involving either mean differences or heteroskedasticity across a pool of n >  groups. we applied statistical theory to determine the specific effect of different means and heteroskedasticity across groups of microarray data on the sign and magnitude of gene-to-gene pearson correlation coefficients obtained from the pool of groups.
 the synthesis of information across microarray studies has been performed by combining statistical results of individual studies, or by combining data from multiple studies into a large pool to be analyzed as a single data set. specific issues relating to data heterogeneity across microarray studies, such as differences within and between labs or differences among experimental conditions, could lead to equivocal results in a melting pot approach.  we applied statistical theory to determine the specific effect of different means and heteroskedasticity across groups of microarray data on the sign and magnitude of gene-to-gene pearson correlation coefficients obtained from the pool of groups. we quantified the biases of the pooled coefficients and compared them to the biases of correlations estimated by an effect-size model. mean differences across the groups were the main factor determining the magnitude and sign of the pooled coefficients, which showed largest values of bias as they approached ± only heteroskedasticity across the pool of groups resulted in less efficient estimations of correlations than did a classical meta-analysis approach of combining correlation coefficients. these results were corroborated by simulation studies involving either mean differences or heteroskedasticity across a pool of n >  groups.  the combination of statistical results is best suited for synthesizing the correlation between expression profiles of a gene pair across several microarray studies. 