BMC Bioinformatics BioMed Central

Methodology article Open Access

A power law global error model for the identification of
differentially expressed genes in microarray data
Norman Pavelka†, Mattia Pelizzola†, Caterina Vizzardelli†,
Monica Capozzoli, Andrea Splendiani, Francesca Granucci and
Paola Ricciardi-Castagnoli*

Address: Department of Biotechnology and Bioscience, University of Milano-Bicocca, Piazza della Scienza 2, 20126 Milan, Italy
Email: Norman Pavelka - [email protected]; Mattia Pelizzola - [email protected];
Caterina Vizzardelli - [email protected]; Monica Capozzoli - [email protected];
Andrea Splendiani - [email protected]; Francesca Granucci - [email protected]; Paola Ricciardi-
Castagnoli* - [email protected]
* Corresponding author †Equal contributors

Published: 17 December 2004 Received: 16 July 2004

Accepted: 17 December 2004
BMC Bioinformatics 2004, 5:203 doi:10.1186/1471-2105-5-203
This article is available from: http://www.biomedcentral.com/1471-2105/5/203
© 2004 Pavelka et al; licensee BioMed Central Ltd.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract
Background: High-density oligonucleotide microarray technology enables the discovery of genes
that are transcriptionally modulated in different biological samples due to physiology, disease or
intervention. Methods for the identification of these so-called "differentially expressed genes"
(DEG) would largely benefit from a deeper knowledge of the intrinsic measurement variability.
Though it is clear that variance of repeated measures is highly dependent on the average expression
level of a given gene, there is still a lack of consensus on how signal reproducibility is linked to signal
intensity. The aim of this study was to empirically model the variance versus mean dependence in
microarray data to improve the performance of existing methods for identifying DEG.
Results: In the present work we used data generated by our lab as well as publicly available data
sets to show that dispersion of repeated measures depends on location of the measures themselves
following a power law. This enables us to construct a power law global error model (PLGEM) that
is applicable to various Affymetrix GeneChip data sets. A new DEG identification method is
therefore proposed, consisting of a statistic designed to make explicit use of model-derived
measurement spread estimates and a resampling-based hypothesis testing algorithm.
Conclusions: The new method provides a control of the false positive rate, a good sensitivity vs.
specificity trade-off and consistent results with varying number of replicates and even using single
samples.

Background problems in microarray data analysis is probably the iden-

DNA microarrays have become common tools for moni- tification of differentially expressed genes (DEG) when
toring genome-wide expression in biological samples har- comparing distinct experimental conditions. In spite of its
vested under different physiological, pathological or biological relevance, there is still no commonly accepted
pharmacological conditions. One of the most challenging way to answer this question.

Page 1 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

An ideal DEG identification method should limit both More recently, Tu et al. [8] empirically modeled the vari-
false positives, i.e. genes wrongly called significant (type 1 ance versus mean dependence from a data set consisting
errors), and false negatives, i.e. genes wrongly called not of ten replicated oligonucleotide microarray experiments.
significant (type 2 errors). To this end, understanding According to the authors, the variance of the genes should
how gene expression values measured in replicated exper- decay exponentially with the mean, but only for moder-
iments are spread around the true expression level of each ately expression values. Taken together, all these aspects
gene, would help to distinguish biologically relevant gene could limit the applicability of these error models.
expression changes from fluctuations due to different
sources of variability that are unrelated to the biological Independently from the choice of the error model,
phenomenon under investigation. Measurement error another point that remains to be faced is on how to esti-
estimates can be obtained in two ways: either by empiri- mate residual error. A discussed by Wright et al. [9] the
cally inferring noise from highly replicated data or by possibilities range between two extremes: either obtaining
deducing noise from a theoretical error model [1]. Espe- a single variance estimate across all genes or obtaining a
cially when the experimental design requires the investi- gene-specific residual variance. In the same paper a hybrid
gation of a high number of conditions, the former strategy approach is proposed in which information from all
is not always feasible, because of the high cost of these genes is used to fit a single linear model from which the
experiments or due to the availability of biological mate- gene-specific variance estimates can be deduced. In the
rial. In addition, there is still a lack of consensus on how present work we chose to follow an approach similar to
gene expression values from replicated experiments the latter.
should be theoretically distributed, which restricts the
application also of the latter strategy. The aim of this study was to use highly replicated micro-
array data to empirically determine the true variance ver-
The most widely used methods for identifying DEG range sus mean dependence that exists in this type of data. This
from purely empirical filtering techniques (e.g. selecting knowledge enabled the proposal of PLGEM as a simple
genes that show a fold change higher than a fixed thresh- but powerful error model. We fitted the proposed model
old) to more sophisticated statistical tests such as the sig- on various data sets without pre-filtering the data, deriv-
nal-to-noise ratio described by Golub et al. [2] or the ing an improved test statistics and identifying DEG even
Significance Analysis of Microarrays (SAM) method by in data sets with very little number of replicates.
Tusher et al. [3]. While empirical filtering techniques rely
on arbitrarily chosen thresholds and are unable to provide Results
any type of control on the significance of the results, the Variance versus mean dependence
more sophisticated statistical tests usually need a high The relationship between measurement variability and
degree of replication in the data to accurately measure average expression values was investigated by means of
gene-specific variability. scatter plots where different measures of spread were dis-
played against different measures of location. For each
In the past years various authors have proposed compet- gene absolute or relative standard deviation was plotted
ing error models for microarray data from which discord- against the mean expression value in either linear or log-
ant implications for the variance versus mean dependence log plots using data from the 16iDC data set (Figure 1).
can be deduced. Chen et al. [4] first proposed a simple Independently from the choice of standard deviation or
Gaussian model, more recently Ideker et al. [5] and Li and inter-quartile range as the estimate of spread and of mean
Wong [6] introduced two-component models containing or median as the location estimate we obtained qualita-
a multiplicative and an additive error term. All of these tively similar plots (data not shown). Log-log plots of
models implicitly or explicitly assume a constant coeffi- both absolute and relative spread estimates revealed a
cient of variation (CV), implying that standard deviation strikingly linear dependency, indicating that measure-
should vary proportionally with the mean. More recently, ment spread could depend on signal location following a
Rocke and Durbin [7] proposed a variation of the two- power law.
component model from which they derived that variance
of repeated microarray measures is a quadratic function of The power law global error model
the mean. Dealing specifically with spotted cDNA micro- Based on the previous observation, we chose to empiri-
array technology Baggerly et al. [1] proposed a beta-bino- cally model measurement noise through linear
mial model, from which it can be derived that variance is regressions:
a second-order polynomial function of the mean. Unfor-
tunately, most of these models are based on theoretical ln( s) = k.ln( x ) + c + ε eqn.1
assumptions that have been verified on simulated data or
on data sets consisting of small numbers of replicates.

Page 2 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

1600
A B8 3
C D2

6 1
1200
2

ln(st.dev.)
4 0
st.dev.

ln(CV)
CV
-2 0 2 4 6 8 10
800
2 -1
1
400
0 -2
-2 0 2 4 6 8 10

0 -2 0 -3
0 1000 2000 3000 4000 5000 6000 7000 0 1000 2000 3000 4000 5000 6000 7000
mean ln(mean) mean ln(mean)

Figure 1 between measurement variability and mean expression level

Relationship
Relationship between measurement variability and mean expression level. For each of the 12488 probe sets dis-
played on Affymetrix MG-U74Av2 chips the standard deviation (st. dev.) or the coefficient of variation (CV = st. dev. / mean) is
plotted against the mean in either linear or log-log plots based on absolute expression values derived from the 16iDC data set.

• Choose a "modeling quantile" q and determine for all

where s and x respectively represent standard deviation
the genes contained in each partition a single "modeling
and mean of repeated measures. Error term ε is the reali-
zation of a random variable E that we will show later to be point" with median of ln( x ) values as the x-coordinate
normally distributed as assumed when fitting a linear and q-th quantile of ln(s) values as the y-coordinate.
model. Inspired by the previous experimental observa-
tions we propose the following power law global error • Finally, find a linear fit through the set of p modeling
model (PLGEM): points using the least-squares method and obtain a slope
k and an intercept c of the resulting regression function.
s( x ) = α ⋅ x β ⋅ eε eqn. 2
Thus, for all possible combinations of p and q a slope kp,q,
Model parameters α and β can be estimated from linear an intercept cp,q and a correlation coefficient r2p,q can be
regression coefficients in 1 in a straightforward way: obtained. Performance of this modeling method was
tested also using different combinations of partitions, in
α = ec eqn. 3 the range between 5 and 500, and quantiles, ranging from
0.01 to 0.99 (Supplementary Table 1). For all 77 analyzed
β=k eqn. 4 combinations of p and q regression lines gave good fit
with the modeling points, with an adjusted r2 that was
PLGEM fitting method always very close to 0.99. In addition, all regression lines
Instead of performing a simple linear fit through the were strikingly parallel as judged by their slopes: 0.81 ±
whole set of points, we preferred to implement a method 0.068 (mean ± sd). The reason for not considering p > 500
that could provide improved model robustness by parti- was that above this number we empirically revealed a
tioning the data to gain local estimates of spread as in poorer modeling quality in terms of correlation coeffi-
Mutch et al. [9]. Most importantly, this method should cients (data not shown), most likely due to the decrease of
also provide the possibility to choose different levels of the number of data points contained in each partition.
confidence when modeling the spread of the data. Note
that Mutch et al. [9] proposed to model within-replicates As a straightforward application of this modeling method,
fold changes as a function of average expression using a PLGEM could be fitted at the 50th-percentile to obtain a
model that was very different from PLGEM. Therefore, the central tendency of standard deviation to be used for
following algorithm was applied: improving test statistics (see next section). Another appli-
cation of this method could be to fit PLGEM at the 5th-
and at the 95th-percentile of standard deviation, to conse-
• Rank genes according to their ln( x ) value and subdivide
quently find the limits of the corresponding 90% empiri-
the overall expression range into a given number p of par-
cal confidence interval of standard deviation.
titions containing an equal number of ranked genes.

Page 3 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

Table 1: Effect of the presence of DEG when applying the classic permutation strategy to the PLGEM-STN statistic.

FPR vs. significance level estimated vs. observed FDR

including DEG excluding DEG including DEG excluding DEG

# of slope intercept adj. R^2 slope intercept adj. R^2 slope intercept adj. R^2 slope intercept adj. R^2
genes in
data set

22300 1.690 3.070 0.856 0.871 -0.114 0.938 0.187 -0.383 0.314 1.090 -0.692 0.826
10000 1.710 2.470 0.881 0.888 -0.083 0.931 0.224 -0.247 0.433 1.040 -0.555 0.824
5000 1.460 1.270 0.880 0.877 -0.147 0.934 0.093 -0.228 0.067 1.080 -0.425 0.836
2500 1.590 1.180 0.888 0.864 -0.155 0.939 0.092 -0.176 0.135 1.110 -0.348 0.853
1500 1.670 0.935 0.908 0.876 -0.166 0.948 0.078 -0.122 0.170 1.130 -0.182 0.880
1000 1.720 0.689 0.874 0.944 0.002 0.956 0.038 -0.122 0.082 0.991 -0.226 0.897
500 1.900 0.263 0.864 0.857 -0.255 0.956 0.062 -0.030 0.307 1.160 0.038 0.909
200 2.490 -0.059 0.875 0.905 -0.378 0.946 0.064 -0.012 0.426 1.050 0.233 0.914

Eight data sets with different percentages of DEG were constructed from the Latin Square data set by keeping the 62 known spiked-in probe sets,
but randomly removing increasing amounts of the remaining probe sets reaching the total indicated in the first column. The null distributions of the
PLGEM-STN statistics were evaluated through the classic permutation strategy either including or excluding DEG. A wide range of significance
levels was used to select DEG and correlation between the FPR and the significance level or between estimated and observed FDR was evaluated
through linear regressions in log-log plots. Table reports slopes, intercepts and adjusted r2 of linear models. See text for details on estimation of
FDR.

In order to verify the feasibility of the former application, cal hypothesis testing, in which we test against the null
fitting of PLGEM on real-life data as well as distribution hypothesis of non-differential expression. First of all, we
properties of the random variable E were investigated by chose to implement as the test statistic the signal-to-noise
analyzing the residuals of the model, i.e. differences ratio (STN) already used by Golub et al. [2], because it
between observed and expected values: explicitly takes unequal variances into account and
because it penalizes genes that have higher variance in
ε g = ln(s g ) − (k10,0.5 ⋅ ln(x g ) + c10,0.5 ) eqn. 5 each class more than those genes that have a high variance
in one class and a low variance in another [11]:
Figure 2A–D shows distribution of residuals εg computed
from the 16iDC data set and its dependency on the rank x (g2) − x (g1)
of mean expression values. Figure 2E–L summarizes STN = eqn. 6
model validation on other two completely unrelated s(g1) + s(g2)
high-density oligonucleotide microarray data sets, the
HG-U133A Leigh syndrome data set (Figure 2E–H) and
x (1) x (2)
the HG-U95Av2 Muscle biopsies data set (Figure 2I–L). where in the original version g and g represent,
For each tested data set an individual model was fitted and respectively, the mean of the replicated expression meas-
a distinct set of parameters α and β was determined. In all
of the three independent data sets measurement variabil- s(1) s(2)
ures of gene g in condition 1 and 2, whereas g and g
ity could be accurately modeled through equation 2, with are the corresponding standard deviations. Instead, we
a power coefficient β that was always between 0 and 1 and propose to use model-derived standard deviation esti-
a random variable E that appeared to be normally distrib- mates predicted by PLGEM in equation 2 for the corre-
uted with zero-mean and constant standard deviation sponding signal mean, rather than data-derived standard
over the whole range of expression values. Of course, deviation values calculated independently from the few
these findings were eventually expected only for q = 0.5, data points that are usually available for every single gene.
and their occurrence demonstrated a goodness of fit of
PLGEM on a series of unrelated real-life data sets. The improvement of the test statistic in ranking DEG was
evaluated as done by Broberg [12] through receiver oper-
Improved test-statistics for detecting differential ator characteristic (ROC) plots on the HG-U133A Latin
expression Square data set, where there is an a priori knowledge on
In order to identify DEG, we implemented the following the truly differentially expressed transcripts. ROC plots
general algorithm derived from the framework of statisti- investigate the relationship between false positive rates

Page 4 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

Figure 2of residuals of PLGEM fitted on three different real-life data sets
Analysis
Analysis of residuals of PLGEM fitted on three different real-life data sets. (A) PLGEM is fitted on the 16iDC data set
(MG-U74Av2) following the method described in the text, setting p = 10 and q = 0.5. (B) Model residuals are plotted as a func-
tion of the rank of the mean absolute expression level. (C) Distribution of pooled residuals. (D) The quantiles of the distribu-
tion of pooled residuals are plotted against the quantiles of a standard normal distribution. The same procedure is applied to
the Leigh syndrome data set (HG-U133A, panels E-H) and to the Muscle biopsies data set (HG-U95Av2, panels I-L).

(FPR) and false negative rates (FNR) at different signifi- statistics for each tested value of n (data not shown). In
cance levels; in this way the performance of the PLGEM- addition, the ROC curve of PLGEM-STN always had the
derived STN statistic (PLGEM-STN) has been compared shortest distance from origin, indicating that it resulted in
with the original STN statistic (CLASSIC-STN) and the the best trade-off between sensitivity and specificity. Inter-
statistic implemented in the commonly accepted Signifi- estingly, improved sensitivity was observed especially
cance Analysis of Microarrays (SAM) DEG identification when the nominal fold change was particularly low (see
method (SAM-STAT). To this purpose Exp01 of the Latin Exp02 vs. Exp01 and Exp14 vs. Exp01).
Square was taken as the baseline to which the remaining
13 experiments were compared. For each comparison Apart from discriminating between significant and not
absolute values of each statistic were ranked in decreasing significant gene expression changes, an optimal test-statis-
order and first n genes selected (where n ranged from 5 to tic should additionally provide an accurate quantification
200). Figure 3 summarizes results only for the most of the actual degree of differential expression. Figure 4
informative comparisons, but in each tested comparison shows that PLGEM-STN outperforms the competing sta-
analysis PLGEM-STN was at least as good as the other two tistics in correlating the value of the statistic with the

Page 5 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

Exp02 vs Exp01 Exp06 vs Exp01

0.003

0.0020
False positive rate

False positive rate

0.002

0.0010
0.001

0.0000
0.000

0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025
False negative rate False negative rate

Exp10 vs Exp01 0.004

0.003
Exp14 vs Exp01
0.0020
False positive rate

False positive rate

0.002
0.0010

0.001
0.0000

0.000

0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025
False negative rate False negative rate

Figure 3
Performance of PLGEM-STN in ranking differentially expressed genes
Performance of PLGEM-STN in ranking differentially expressed genes. ROC plots were used to compare the sensi-
tivity vs. specificity trade-off of the following three statistics: PLGEM-STN (black), CLASSIC-STN (blue) and SAM-STAT (red).
SAM-STAT values were obtained using the R package "siggenes" [21]. Absolute values of the corresponding statistics were
sorted in decreasing order, first n genes were selected (where n ranged from 5 to 200) and false positive and false negative
rates were evaluated on the HG-U133A Latin Square dataset. Note that, while the transcripts in Exp02 (Exp14) are spiked-in
at twice (half) the concentration than in Exp01, in both Exp06 vs. Exp01 and in Exp10 vs. Exp01 comparisons the nominal fold-
change of spiked-in transcripts ranged from 32 to 512.

nominal concentration variation of the known Latin

Square DEG; this was particularly true for the most
extreme variations.

Page 6 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

Nevertheless, we can foresee that the classic permutation

PLGEM STN
40

40
strategy may not be optimal for estimating the actual FPR
0

0
when the test-statistic makes use of a global error model
such as PLGEM. We can in fact hypothesize that measure-
-40

-40
ment spread of DEG may not be accurately described by
-10 -5 0 5 10 0.1 1.0 10.0 500.0 means of a global error model that was designed to
describe signal variability in absence of differential expres-
100
CLASSIC STN

sion. To test this hypothesis we compared the correlation

50
between the expected significance level and the observed
0

FPR using PLGEM-STN and the classic permutation strat-

-50
-100

egy either including or excluding DEG during the permu-

-10 -5 0 5 10 0.1 1.0 10.0 500.0 tation step. To this end, data sets containing different
percentages of DEG were obtained by merging the 62
150
SAM's statistic

known DEG of the Latin Square data set with differently

40

sized random samples of not DEG extracted from the

0

same data set. As predicted, the presence of DEG during

-40

the permutation step caused the significance level to be

-150

less correlated with the observed FPR and this correlation

-10 -5 0 5 10 0.1 1.0 10.0 500.0
worsened with increasing percentages of DEG (Table 1).
Nominal Log Ratio Concentration (pM) This lack of correlation was dramatically amplified when
expected and observed numbers of false positives were
statistics
Correlation
inal
Figure
concentration
4 between variation
the value
in comparison
of PLGEM-STNto competing
and the nom-
Correlation between the value of PLGEM-STN and divided by the number of selected genes to obtain an over-
the nominal concentration variation in comparison simplified estimate of the false discovery rate (FDR) and
to competing statistics. Exp01 of the Latin Square data the observed FDR. Conversely, when DEG were omitted
set was taken as the baseline to which the remaining 13 during the permutation step the correlation between esti-
experiments were compared. The observed values of the mated and observed FPR or estimated and observed FDR
indicated statistics are plotted against the nominal log ratio, was sensibly higher for each tested percentage of DEG. We
deduced from the known spiked-in concentrations (left pan- hereby by no means claim that this FDR estimate is the
els). A nominal log ratio of 0 is assumed for the remaining most accurate. A more appropriate relationship between
transcripts and a box-plot of their corresponding values of FPR and FDR can be found in the paper by Storey and Tib-
the indicated statistics is superimposed to the plot. For those shirani [13]. Nevertheless, the explicit control of the FDR
cases where one of the two known spiked-in concentrations
goes beyond the scope of the present paper.
is 0, the value of the statistic is instead plotted against the
non-null concentration (right panels). Red and black dots
represent transcripts that are present in Exp01 or in the Since in real-life data sets true DEG are unknown in
remaining 13 experiments, respectively. advance, we propose the following resampling-based
method to obtain the null distribution of not DEG when
comparing n1 replicates of condition A with n2 replicates
of condition B:
Identification of differentially expressed genes
A resampling-based method for estimating the null distribution • Artificial condition A* is obtained by randomly sam-
Though ranking of genes based on the absolute value of pling with replacement n1 indexes corresponding to the
their test-statistic has been proven to be an effective replicates of only one experimental condition. If availa-
method for selecting DEG, an even more useful way ble, chose the condition with the highest number of
would be to compare the observed statistic with its null replicates;
distribution (the distribution of values of the statistic that
are expected by chance for a not differentially expressed • Similarly sample n2 values from the same set to obtain
gene), in order to control the FPR. indexes of artificial condition B*;

A classic approach to empirically obtain the null distribu- • Compute resampled test-statistics between A* and B* at
tion of a test-statistic is running a series of random permu- each cycle.
tations of the chip indexes of the full data set and re-
computing the test-statistics at each permutation. Permu- The previous resampling should be repeated a sufficiently
tated test-statistics can then be pooled and significance large number of times – as large as possible compared to
thresholds (i.e. expected false positive rates) are found as the total number of possible combinations and compati-
specific quantiles of the null distribution. bly with available computational resources – and the

Page 7 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

pled PLGEM-STN values have a good concordance with

Permutated PLGEM-STN (mean quantiles)

the mean quantiles of the classically permutated statistics

averaged over the 13 comparisons, implying that no dif-
4

ferences are expected also in the gene selection step.

In accordance with the previous observations, the ROC

2

curve of the resampling method applied to the PLGEM-

STN statistic was not significantly different from the ROC
curve of the classic permutation strategy (excluding DEG)
0

applied to the same statistic on the Latin Square data set

(data not shown). Conversely, ROC curves of the classic
permutation strategy (including DEG) applied to the
-2

CLASSIC-STN statistic and of the SAM method gave

poorer performance similarly to the results in Figure 3
(data not shown).
-4

1.0
-3 -2 -1 0 1 2 3

0.8
Resampled PLGEM-STN (quantiles)
Figure
Comparison
tion of the
5 PLGEM-STN
of two methods
statistic
for inferring the null distribu- 0.6
Frequency

Comparison of two methods for inferring the null dis-

tribution of the PLGEM-STN statistic. The classic per-
mutation strategy (excluding DEG) was performed for each
0.4

comparison in the Latin Square data set and the quantiles of

the distribution of PLGEM-STN values were averaged over
the 13 comparisons. The mean quantiles of the permutated
0.2

statistics are plotted against the quantiles of the distribution

of PLGEM-STN values obtained through the proposed resa-
mpling approach performed on the same data set but includ-
0.0

ing DEG.

2 4 6 8 10

Counts

resampled test-statistics finally pooled. In our opinion Figure are

Consistency
thereof 6 analyzed
of findings when different replicates or numbers
resampling the expression values from only one experi- Consistency of findings when different replicates or
mental condition, rather than permutating indexes of numbers thereof are analyzed. Ten reduced data sets
both conditions, makes more sense with this particular were constructed by removing all possible combinations of 1
or 2 replicates of LPS-stimulated DC from the data set. The
statistic, because in this way we avoid merging true and
plot shows the cumulative frequency at which the probe sets
false null hypothesis. Note that when more than one con- are consistently selected in the 16iDC+LPS and in the ten
dition (all with the same number of replicates) are to be reduced data sets by the following methods: the resampling
compared to a common baseline, the distribution of resa- approach applied to PLGEM-STN (filled squares), the permu-
mpled test-statistics needs to be determined only once, tation strategy applied to CLASSIC-STN (filled circles) and
obviously providing a computational advantage. As a test the SAM method (open triangles). In case of PLGEM, a single
of substantial equivalence between this resampling model was fitted on the common 16iDC baseline data set.
method and the classic permutation strategy (excluding The cumulative frequencies are normalized with respect to
DEG), we compared the distribution of the permutated the total number of probe sets identified by the correspond-
and of the resampled PLGEM-STN test-statistics in Q-Q ing method.
plots. The distribution of the PLGEM-STN resampled
from Exp01 of the Latin Square data set was almost iden-
tical with the distributions of permutated PLGEM-STN
obtained with the classic strategy from each comparison Increased robustness to varying number of replicates
with the remaining 13 experimental conditions (data not Another appealing feature of an optimal DEG identifica-
shown). Figure 5 shows that the quantiles of the resam- tion method is that it should provide consistent results

Page 8 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

when different replicates of a same data set or different Discussion

numbers thereof are analyzed. We therefore compared the PLGEM accurately describes GeneChip data variability
performance of our resampling approach applied to the In the present work we described a new global error
PLGEM-STN statistic (method 1) with SAM (method 2) model for microarray gene expression data that describes
and with the classic permutation strategy applied to the measurement variability with the same degree of accuracy
CLASSIC-STN statistic (method 3). The number of availa- over the whole dynamic range of values and that can be
ble replicates for each experimental condition in the Latin fitted at any desired quantile of spread. PLGEM has
Square data set was unfortunately too small to investigate proven to correctly model signal standard deviation, in
this particular task. We therefore took advantage of the spite of the presence of different sources of variability, e.g.
16iDC+LPS data set, where the first sixteen columns can biological variability as well as the use of different target
be considered as the baseline condition for the remaining preparation protocols or of different chips. Moreover,
four experimental replicates. We then constructed a series PLGEM has shown to be able to deal with the great varia-
of reduced data sets in which the baseline columns were bility that exists at low expression levels while at the same
kept constant while all possible combinations of 1, 2 or 3 time considering the significant relative reproducibility of
replicates of LPS-stimulated DC were systematically highly expressed genes. Previously proposed error models
deleted from the 16iDC+LPS data set, reaching a total of assumed that measurement spread depended on signal
fifteen distinct data sets including the original one. Since location following different mathematical relationships,
methods 2 and 3 are not applicable on the four reduced but none of them was based on a power law thus far. Anal-
data sets containing single samples for the LPS experimen- ysis of the residuals showed a good fit of PLGEM to a
tal condition, only the eleven data sets with at least two number of high-density oligonucleotide microarray data
replicates were used for comparison purposes. Since the sets, with model parameters being very similar to each
sixteen baseline columns are identical in each reduced other even when dealing with RNA samples coming from
data set, PLGEM parameters were determined only once completely different biological sources and analyzed on
on this common baseline condition. Significance levels different array layouts. This suggests that PLGEM could
used by each method in all eleven data sets were empiri- represent a general Affymetrix GeneChip measurement
cally selected in order to achieve a similar number of sig- noise model. Even though scaled MAS5 Signals gave satis-
nificant genes (ca. 500 probe sets) in the full data set, i.e. factory modeling results, a further improvement could be
the one containing all available replicates. Thus, for each achieved by using other emerging gene expression indices
method eleven lists of identified DEG were obtained and [6,14] or more sophisticated normalization techniques,
the consistency between these lists was evaluated by e.g. quantile normalization [15]. Interestingly, if the same
counting the number of times each probe set was selected, evaluation of sensitivity vs. specificity using ROC plots on
giving a probes set count between 1 and 11. In Figure 6 we the Latin Square data set was done using GCRMA expres-
compared the three distinct cumulative frequency curves sion values [16], the results were even more striking than
for each method, which show the percentage of identified using MAS5 Signals (data not shown). Further studies will
DEG that were selected at least a given number times. be needed to assess if PLGEM is also able to deal with data
While method 2 and 3 gave similar results, the method coming from microarray technologies others than
proposed in the present work identified a larger number Affymetrix GeneChips.
of probe sets in a larger number of lists.
Interestingly, model parameter β was found to be quite
We finally evaluated the possibility of applying our stable and comprised between 0 and 1 in all analyzed data
method also to data sets where one of the experimental sets. It is noteworthy that for β ∈ (0:1) absolute variability
conditions was investigated only with a single sample increases with growing expression values, while relative
without replication. To this end, we used the remaining variability decreases (compare panel B with panel D of
four reduced data sets that could not be used in the previ- Figure 1). On the other hand, none of the models men-
ous comparison. In this case, the same PLGEM parameters tioned in the background section seem to agree with these
derived from the sixteen baseline columns were applied to experimental observations. Formal statistical reasoning
each of the single LPS-treated DC sample to obtain an esti- could unravel the underlying theoretical error model that
mate of standard deviation associated to each gene expres- leads to the power law relationship that was observed to
sion value, treated here as if it was a mean value from a be at the basis of the variance versus mean dependence in
larger group of values. Interestingly, when results replicated microarray data.
obtained through this procedure were compared to the
previously described results a comparable number of DEG A PLGEM-based method successfully detects differential
was identified and only one probe set was newly detected expression
in comparison to the previously identified ones (data not In spite of the lack of a theoretical statistical model, the
shown), arguing for a good consistency of results. empirical model presented here has proven its

Page 9 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

applicability in the identification of DEG, providing Conclusions

improved results under a wide range of different testing The proposed DEG identification method provides a
conditions. In comparison to other commonly used DEG direct control of the FPR and an indirect control of the
identification methods, the proposed approach demon- FDR. Moreover, as tested on the Latin Square data set, our
strated improved specificity and sensitivity on the Latin method improved the specificity vs. sensitivity trade-off in
Square data set and robustness to decreasing number of comparison to other commonly applied DEG selection
replicates on the 16iDC+LPS data set. The good perform- techniques. It finally showed an increased robustness
ance of our proposed method is reasonably due to the fact when different replicates or numbers thereof are analyzed,
that it relies on a global error model. As an example, when giving consistent results even in data sets containing sin-
the classic permutation strategy is applied to the CLASSIC- gle samples. In conclusion, the global error model pre-
STN statistic or when the SAM method is used, the sented here may facilitate the analysis of microarray gene
selected genes are apparently more dependent on the expression data by discriminating information from
number and identity of the replicates than when our pro- noise, and thus possibly helping the formulation of new
posed approach is used. We hypothesize that, when no hypothesis concerning gene functions.
error model is assumed and a small number of replicates
is present in the data set, the probability of observing for Methods
some genes coincidently very similar (or very dissimilar) Data sets
values increases, thus leading to an underestimation (or 16iDC
overestimation) of the standard deviation and a conse- RNA was harvested from ten biological samples of
quent overestimation (or underestimation) of the test sta- unstimulated immature mouse dendritic cells (DC), each
tistic, finally leading to false positives (or false negatives). extracted from an independent batch of cells. One opera-
tor prepared the biotin-labeled cRNA for hybridization
Interestingly, when the performance of our method was from three of the ten RNA samples, a second operator pre-
compared on a data set of DC stimulated for 24 hours pared the remaining seven. While operator 1 applied the
with LPS, SAM showed a decreased sensitivity in identify- total RNA protocol to all of its three samples, operator 2
ing down-regulated genes when the number of LPS repli- applied the purified mRNA protocol to five of its seven
cates was low (data not shown). Under these samples and the total RNA protocol to the remaining two.
experimental conditions DC undergo a process known as Two of the three cRNA samples prepared by operator 1
maturation, which is a specialized form of cellular and four of the seven cRNA samples prepared by operator
differentiation, for which both up- and down-regulation 2 have been hybridized twice; therefore, a total of 16 MG-
of gene expression is expected [17,18]. We speculate that U74Av2 GeneChips (Affymetrix, Santa Clara, CA) have
SAM did not select these genes, because of the combina- been employed.
tion of two effects. First of all, down-regulated genes are
expected to have lower and therefore intrinsically more Leigh syndrome
variable expression values in the four LPS replicates than Eight RNA samples were harvested from human fibroblast
in the sixteen replicates of immature DC. When, in addi- cell lines each deriving from a distinct Leigh syndrome
tion, the number of LPS replicates becomes too low, SAM patient [19,20] and individually hybridized on HG-
filters these genes out to control the FDR. In agreement U133A GeneChips (Affymetrix).
with this hypothesis SAM was perfectly able to identify
down-regulation when the full data set was used (data not Muscle biopsies
shown). Four individual and two pooled RNA samples from
human muscle biopsies of sixteen healthy young male
The gene selection method proposed in the present work donors were hybridized on six HG-U95Av2 GeneChips
does not provide a direct control on the FDR, but the sig- (Affymetrix). This data set was downloaded from [21],
nificance level has been proven to be a direct estimate of experiment code: GSE80 [22].
the FPR. Thus, if a significance level of 0.001 is used and
12488 probe sets are displayed on the MG-U74Av2 chip, Latin Square
12–13 genes are expected to be selected by chance in cases This data set consists of 3 technical replicates of 14 sepa-
where all genes are in fact not differentially expressed. rate hybridizations (named Exp01–14) of 42 spiked tran-
Therefore, a researcher can test how many genes would be scripts in a complex human background at concentrations
selected over a range of different significance levels and ranging from 0.125 pM to 512 pM. Thirty of the spikes are
chose the one that results in the most acceptable compro- isolated from a human cell line, four spikes are bacterial
mise between number of selected genes and estimated controls, and eight spikes are artificially engineered
FPR. sequences believed to be unique in the human genome.
Further details on the design of the Latin Square data set

Page 10 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

can be found at [23]. Considering the redundancy of Performance of modeling method using different combinations of
some probe sets, there are a total of 62 distinct probe sets parameters p and q. The modeling method described in this study was
tested on the 16iDC data set using different combinations of partitions (5,
designed to match the 42 spiked transcripts.
10, 20, 50, 100, 200 and 500), and quantiles (0.01, 0.02, 0.05, 0.1,
0.2, 0.5, 0.8, 0.9, 0.95, 0.98 and 0.99). For all 77 analyzed combina-
16iDC+LPS tions of p and q regression lines were fitted to the data as described in the
This data set consists of the same samples of the 16iDC text. Goodness of fit was evaluated from the resulting slope (panel A),
data set, but includes additional four samples as a second intercept (panel B) and adjusted r2 (panel C).
experimental condition. To this end dendritic cells were Click here for file
[http://www.biomedcentral.com/content/supplementary/1471-
stimulated to mature with lipopolysaccharide (LPS) for 2105-5-203-S1.xls]
24 hours. Two independent biological samples were har-
vested and individually processed by the same two opera-
tors that prepared the samples for the 16iDC data set: one
applied the total RNA protocol, the other one applied the Acknowledgements
purified mRNA protocol. Each cRNA sample was This research was supported by AIRC (Italian Association for Cancer
hybridized twice, thus using a total of four Affymetrix MG- Research) and by a grant from the CARIPLO Foundation ("Development of
U74Av2 chips. Functional Genomics and Bioinformatics platforms aimed to foster novel
approaches in immunotherapy and molecolar diagnostics"). Authors would
Software like to kindly acknowledge Stefano Monti (Broad Institute, Cambridge, MA,
USA) and Peter J. Park (Harvard Medical School, Boston, MA, USA) for
All chips mentioned in the present study were hybridized
helpful discussion and Ottavio Beretta, Gianpiero Cattaneo, Maria Foti,
and scanned following Affymetrix recommendations and Giorgio Moro and Ettore Virzi (University of Milano-Bicocca, Milan, Italy)
MicroArray Suite 5.0 (MAS5) was used as the image acqui- for critical reading of the manuscript. We are also grateful to Valeria
sition and analysis software. All data sets used passed Tiranti, Rossana Mineri and Massimo Zeviani (Unit of Molecular Neuroge-
quality control tests and probe set signals were scaled so netics, National Neurological Institute "Carlo Besta", Milano, Italy) for
that the 4%-trimmed mean of all expression values of kindly providing the HG-U133A Leigh syndrome data set.
each chip was equal to a predefined reference intensity
(called TGT) following manufacturer's recommendations: References
1. Baggerly KA, Coombes KR, Hess KR, Stivers DN, Abruzzo LV, Zhang
W: Identifying differentially expressed genes in cDNA micro-
TGT = 100 for MG-U74Av2 and HG-U133A chips and array experiments. J Comput Biol 2001, 8:639-659.
TGT = 500 for HG-U95Av2 chips. 2. Golub TR, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov
JP, Coller H, Loh ML, Downing JR, Caligiuri MA, Bloomfield CD,
Lander ES: Molecular classification of cancer: class discovery
All procedures for fitting PLGEM, for calculating observed and class prediction by gene expression monitoring. Science
PLGEM-based signal-to-noise ratios (STN), for obtaining 1999, 286:531-537.
3. Tusher VG, Tibshirani R, Chu G: Significance analysis of micro-
expected PLGEM-STN through the resampling-based arrays applied to the ionizing radiation response. Proc Natl
approach and for comparing observed with expected STN Acad Sci U S A 2001, 98:5116-5121.
values have been implemented as R functions [24] and 4. Chen Y, Dougherty ER, Bittner ML: Ratio based decisions and the
quantitative analysis of cDNA microarray images. J Biomed
will be soon submitted for integration into the Biocon- Opt 1997, 2:364-374.
ductor project [25]. 5. Ideker T, Thorsson V, Siegel AF, Hood LE: Testing for differen-
tially-expressed genes by maximum-likelihood analysis of
microarray data. J Comput Biol 2000, 7:805-817.
Authors' contributions 6. Li C, Wong WH: Model-based analysis of oligonucleotide
NP conceived the study and drafted the manuscript. MP arrays: expression index computation and outlier detection.
Proc Natl Acad Sci U S A 2001, 98:31-36.
wrote the software, participated in the design of the study 7. Rocke DM, Durbin B: A model for measurement error for gene
and in the editing of the manuscript. CV performed the expression arrays. J Comput Biol 2001, 8:557-569.
8. Tu Y, Stolovitzky G, Klein U: Quantitative noise analysis for
microarray experiments, participated in the design of the gene expression microarray experiments. Proc Natl Acad Sci U
study and the editing of the manuscript. MC participated S A 2002, 99:14031-14036.
in the microarray experiments, AS participated in the 9. Wright GW, Simon RM: A random variance model for detec-
tion of differential gene expression in small microarray
design of the algorithms, FG and PRC coordinated the experiments. Bioinformatics 2003, 19:2448-2455.
study. All authors read and approved the final manuscript. 10. Mutch DM, Berger A, Mansourian R, Rytz A, Roberts MA: The limit
fold change model: A practical approach for selecting differ-
entially expressed genes from microarray data. BMC
Additional material Bioinformatics 2002, 3:17.
11. Cancer Genomics Software: GeneCluster 2.0b Reference
Manual [http://www.broad.mit.edu/cancer/software/genecluster2/
Additional File 1 gc_ref.html]
12. Broberg P: Statistical methods for ranking differentially
expressed genes. Genome Biol 2003, 4:. epub May 29
13. Storey JD, Tibshirani R: Statistical significance for genomewide
studies. Proc Natl Acad Sci U S A 2003, 100:9440-9445.
14. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ,
Scherf U, Speed TP: Exploration, normalization, and summa-

Page 11 of 12
(page number not for citation purposes)
BMC Bioinformatics 2004, 5:203 http://www.biomedcentral.com/1471-2105/5/203

ries of high density oligonucleotide array probe level data.

Biostatistics 2003, 4:249-264.
15. Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of
normalization methods for high density oligonucleotide
array data based on bias and variance. Bioinformatics 2003,
19:185-193.
16. Wu Z, Irizarry RA, Gentleman R, Murillo FM, Spencer F: A Model
Based Background Adjustment for Oligonucleotide Expres-
sion Arrays. Johns Hopkins University, Dept. of Biostatistics Working
Papers, Working Paper 1 [http://www.bepress.com/jhubiostat/paper1].
(May 28, 2004)
17. Granucci F, Vizzardelli C, Virzi E, Rescigno M, Ricciardi-Castagnoli P:
Transcriptional reprogramming of dendritic cells by differ-
entiation stimuli. Eur J Immunol 2001, 31:2539-2546.
18. Granucci F, Vizzardelli C, Pavelka N, Feau S, Persico M, Virzi E,
Rescigno M, Moro G, Ricciardi-Castagnoli P: Inducible IL-2 pro-
duction by dendritic cells revealed by global gene expression
analysis. Nat Immunol 2001, 2:882-888.
19. Tiranti V, Hoertnagel K, Carrozzo R, Galimberti C, Munaro M, Gra-
natiero M, Zelante L, Gasparini P, Marzella R, Rocchi M, Bayona-Bafa-
luy MP, Enriquez JA, Uziel G, Bertini E, Dionisi-Vici C, Franco B,
Meitinger T, Zeviani M: Mutations of SURF-1 in Leigh disease
associated with cytochrome c oxidase deficiency. Am J Hum
Genet 1998, 63:1609-1621.
20. Tiranti V, Jaksch M, Hofmann S, Galimberti C, Hoertnagel K, Lulli L,
Freisinger P, Bindoff L, Gerbitz KD, Comi GP, Uziel G, Zeviani M,
Meitinger T: Loss-of-function mutations of SURF-1 are specif-
ically associated with Leigh syndrome with cytochrome c
oxidase deficiency. Ann Neurol 1999, 46:161-166.
21. Gene Expression Omnibus database [http://
www.ncbi.nlm.nih.gov/geo]
22. Welle S, Brooks AI, Thornton CA: Computational method for
reducing variance with Affymetrix microarrays. BMC Bioinfor-
matics 2002, 3:23.
23. Affymetrix – Latin Square Data [http://www.affymetrix.com/
support/technical/sample_data/datasets.affx]
24. The R Project for Statistical Computing [http://www.r-
project.org]
25. Bioconductor [http://www.bioconductor.org]

Publish with Bio Med Central and every

scientist can read your work free of charge
"BioMed Central will be the most significant development for
disseminating the results of biomedical researc h in our lifetime."
Sir Paul Nurse, Cancer Research UK

Your research papers will be:

available free of charge to the entire biomedical community
peer reviewed and published immediately upon acceptance
cited in PubMed and archived on PubMed Central
yours — you keep the copyright

Submit your manuscript here: BioMedcentral

http://www.biomedcentral.com/info/publishing_adv.asp

Page 12 of 12
(page number not for citation purposes)