September 4, 2009

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

1

QUALITY CONTROL AND REPRODUCIBILITY IN DNA
MICROARRAY EXPERIMENTS
ANDRÉ FUJITA1
andrefujita@riken.jp
MARIA C. GALVÃO4
mchrisgalvao@gmail.com

JOÃO R. SATO2
joao.sato@ufabc.edu.br
MARI C. SOGAYAR 3
mcsoga@iq.usp.br

FERNANDO H.L. DA SILVA3
lojudice@iq.usp.br
SATORU MIYANO1,5
miyano@ims.u-tokyo.ac.jp

1 Computational Science Research Program, RIKEN, 2-1 Hirosawa, Wako, Saitama,351-0198, Japan
2 Center of Mathematics, Computation and Cognition, Universidade Federal do ABC,

Rua Satna Adélia, 166 - Santo André
3 Chemistry Institute, University of São Paulo, Av. Lineu Prestes, 748 - São Paulo,

05508-000, Brazil
4 Orthodontic Department, University Methodist of São Paulo, Rua Alfeu Tavares,

149 - São Bernardo do Campo, 09641-000, Brazil
5 Human Genome Center, Institute of Medical Science, University of Tokyo, 4-6-1

Shirokanedai, Minato-ku, Tokyo 108-8639, Japan
Biological experiments are usually set up in technical replicates (duplicates or triplicates) in order to ensure reproducibility and, to assess any signiﬁcant error introduced
during the experimental process. The ﬁrst step in biological data analysis is to check thetechnical replicates and to conﬁrm that the error of measure is small enough to be of noconcern. However, little attention has been paid to this part of analysis. Here, we propose a general process to estimate the error of measure and consequently, to provide aninterpretable and objective way to ensure the technical replicates’ quality. Particularly,
we illustrate our application in a DNA microarray dataset set up in technical duplicates.
Keywords: reproducibility; quality control; technical replicate; Dahlberg’s error; DNA
microarrays

1. Introduction
Technical replicates are mandatory in most biosciences experiments in order toensure the results consistency and to identify and avoid possible measurement errors
derived from methodological/technical processes. However, little attention has been
paid to the analysis and quality control among replicates (for convenience, along
all the text, we will refer to “technical replicate” simply as “replicate”). In general,
discrimination between acceptable and poor replicates is carried out in a subjective
fashion, without interpretable or well-deﬁned criteria and the replicates quality
control is not judiciously described, with little attention being paid to this initial
step of data analysis.
It is a common sense among biologists to use a pre-deﬁned value to distinguishacceptable replicates from poor replicates. For example, performing a Real Time

September 4, 2009

2

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

A. Fujita et al.

RT-PCR experiment in duplicates, one may assume that a diﬀerence between thereplicates may not be greater than a given a priori set ∆. Despite that, this predeﬁned number ∆ does not diﬀerentiate between highly and lowly expressed genes,
i.e., it is not proportional to the gene expression value. To illustrate this limitation,
suppose ∆ = 0.5 and two genes, one with expression values equal to, respectively,
10 and 30. The proportion between 0.5 in 10 and 0.5 in 30 is totally diﬀerent, i.e.,∆ may be more restrictive for low expression genes than high expression ones. In
addition, it is diﬃcult to statistically quantify the error of this pre-deﬁned threshold,which varies from one operator to the other.
A second problem is the analysis of the microarray data. It has been describedthat some DNA microarrays such as Aﬀymetrix, Agilent and Codelink [1, 3] provide
Pearson correlation coeﬃcients between replicates greater than 0.9. For other platforms, such as cDNA microarrays or the Mergen platform, the Pearson correlationcoeﬃcient between technical replicates varies from a low value of 0.5 and a highvalue of 0.95 [3, 12, 13]. For a detailed review about reproducibility in microarrays,
see Draghici et al. [7].
As described above, microarray reproducibility is usually measured by applying
Pearson’s correlation [16], in which a value close to one indicates good reproducibility, otherwise, a bad reproducibility. However, Pearson’s correlation assumes thatthe variance along the data is equal, when it is known that the variance along the
spots in a microarray varies, i.e., there is heteroscedasticity [2, 26]. Moreover, Pearson’s correlation is a measure of proportionality, and not a measure of how muchthe values of the spots from the second microarray are similar to those of the ﬁrst
one (reproducibility). Mean and amplitude biases (systematic errors) cannot be detected by correlation coeﬃcient, therefore, Pearson’s correlation cannot be applied
to verify reproducibility. Kim et al. [15] have proposed the use of Spearman’s correlation, nevertheless, analogously to Pearson’s, the former is a measure of associationand not of error.
In addition, correlation’s measures do not identify the poor quality spots or the
reproducibility of speciﬁc spots, but only provide the general association between
the replicated microarrays.
Here, we suggest a general method to estimate the measurement error basedon the concept of Dahlberg’s error [5], which may be applied in most biologicalexperiments involving technical replicates. Moreover, we present a solution which
overcomes some well-known limitations of Dahlberg’s error, such as homoscedasticity and incorporation of bias (systematic errors). We also apply the proposed
method to actual DNA microarray data in order to illustrate usefulness of this approach. To this end, we model the heteroscedasticity, by providing a quality control
criterion for each spot, based on the replicated data.

September 4, 2009

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

Quality Control in Biological Experiments

3

2. Materials and Methods
Primarily, we will describe an error of measure’s estimator based on the Dahlberg’s
method and its limitations. We then introduce an estimation, based on Support
Vector Regression, in order to overcome these limitations, and an algorithm to deal
with the heteroscedasticity problem in DNA microarray data.
2.1. Dahlberg’s Error (D.E.)
Consider the following model:
Zij = µi + εij(1)where Zij is the measure obtained in one biological experiment, i is the sample
index i = 1, . . . , N , j is the replicate number (j = 1, 2 in the case of duplicates), µi
is the unknown true value of the measure and εij is the error of measure.
For the error of measure, assume that E(εij ) = 0 and V ar(εij ) = δε2 . Thus, one
possible quantiﬁcation of the measure’s quality is the standard deviation of εij , i.e.,δε . In other words, the lower is the standard deviation of the error of measure, themore reproducible the method is.
Considerdi = Zi2 − Zi1(2)

V ar(di ) = V ar(εi2 − εi1 ) = 2δε2(3)Therefore,

Assuming that there is no bias (systematic error), one intuitive estimator for2δε2 is2δ̂ε2 =

N
∑
d2
i

N
i=1

(4)

√
∑N d2iThe quantity δ̂ε =i=1 2N is exactly the Dahlberg’s formula proposed in1940 in order to estimate the error of measure in cephalometric studies [5, 10, 11].
This estimator for the standard deviation of the error of measure is widely used in
Orthodontics [5, 6, 10, 11, 14, 20, 22] and may be interpreted as root of squared
error’s average.
One may use this standard deviation of the error of measure in order to check
whether the replicates are similar enough or not. It is known that approximately
95% of the values of a random sample generated from a normal distribution has
a mean between µ − 1.96δε2 and µ + 1.96δε2 . Consequently, one criterion to verify
whether replicates are similar or not may be:September 4, 2009

4

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

A. Fujita et al.

T =1 1.96δε1.96δε
(+)>α
2 |Zi1 ||Zi2 |(5)This quantity T indicates that in approximately 95% of the sample the ratio
between the error of measure and the observed measure is lower than T , in which
T is the proportion of error related to the measured value, allowing evaluation of
the performed error.
If T is greater than a deﬁned α, the replicate may be excluded and the experimentshould be repeated. Notice that T is proportional to the measured data.
Unfortunately, the Dahlberg’s error is extremely aﬀected by systematic errors. Notice that any bias between the two measures is incorporated. Moreover,Dahlberg’s error assumes equal means and variances between both measures. Therefore, Dahlberg’s error does not discriminate between systematic error (biases inmeasurement which leads to measured values systematically higher or lower than
the true value) and random error (unpredictable ﬂuctuations in the measurements),
rendering the interpretation of the results very diﬃcult.
In order to overcome these limitations, we suggest an approach based on Support
Vector Regression.
2.2. Support Vector Regression (SVR)
SVR is a robust regression developed by Vapnik and Lerner [23] and Vapnik andChervonenkis [24] and recently applied to Bioinformatics [8, 9, 17].
Let {(x1 , y1 ), . . . , (xi , yi )} ⊂ R × R be the values obtained from biological experiments performed in duplicates, where y is the replicate of x.
In ε-insensitive SVR [25], it is estimated by a function f (x) that has at the mostε deviation from the yi for all the data, and is as ﬂat as possible.
In other words, the intuitive idea is to deﬁne a tube of radius ε around theregression, where ε > 0, and no error is computed if y lies inside the tube. Therefore,
outliers are naturally excluded from the regression computation (see Figure 1).
More technically, in the case of linear functions f :f (x) = hwt xi + b
with w ∈ R , b ∈ R.
Flatness in (6) means small w, i.e., minimize ||w||2 .
Minimize12
kwk
2
Constrained to

(6)

n

{

yi − (wt xi ) − b ≤ ε(wt xi ) + b − yi ≤ ε

(7)

(8)

September 4, 2009

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

Quality Control in Biological Experiments

Fig. 1.

5

An illustrative scheme of the linear SVR.

Notice that in (8) there is a function f which, with ε precision, approximates all
pairs (xi , yi ). In the cases where it is necessary to allow for some errors, the problem
can be reformulated to [25]:Minimize∑1
2
kwk + C(ξi + ξi∗ )(9)2Constrained tot yi − (w xi ) − b ≤ εit(w xi ) + b − yi ≤ ε + ξi∗ξi , ξi∗ ≥ 0(10)where ξi , ξi∗ are slack variables, the constant C > 0 is the trade-oﬀ between theamount up to which deviations larger than ε are tolerated, maintaining the ﬂatnessof f .
For more details about numerical computation and theoretical SVR information,see Smola and Schölkopf [21].
2.3. Modeling DNA microarray data
As described in the Introduction, for microarray data, it is known that the variancevaries along the spots due to technical problems such as hybridization eﬃciency,
probe sequence, background ﬂuorescence, signal quantiﬁcation procedures [27, 28],
therefore, application of the Dahlberg’s formula is not straightforward. In order toovercome this problem, we suggest the following algorithm:
Let X and Y be two DNA microarrays, with Y being the replicate of X.

September 4, 2009

6

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

A. Fujita et al.

(1) Perform a non-linear regression which is robust to outliers, namely Support
Vector Regression [8] between log(X) and log(Y ), i.e., log(Y ) = f (log(X)) + ε1 .
Notice that the logarithm was calculated due to the high variance observed in
microarray data. This is a common practice in microarray data analysis. For
other biological replicates, which do not present high variance, such as RealTime RT-PCR, it is not necessary to apply logarithm.
(2) Apply again the Support Vector Regression between ε21 and log(X), i.e., ε̂21 =f (log(X)) + ε2 .
(3) Calculate δ̂ = f (log(X)), which is exactly the error of measure. Notice that with2
this process, we obtain one δ̂i for each spot i = 1, . . . , N , where N is the numberof spots in the microarray.
1.96eδi
1.96eδi(4) Calculate Ti = f (log(X+ f (log(X, i = 1, . . . , N , where N is the numberi ))Xii ))Yiof spots in the microarray.
It is important to note that this method is based on the normality of the residues,
therefore, this condition must be checked.
This SVR-based method may be applied to any replicated data such as RealTime RT-PCR, DNA microarrays, protein quantiﬁcations etc.
2.4. DNA Microarray
2.4.1. Cell Lysis and RNA Extraction
Cell cultures were lysed and their RNA extracted using the Illustra RNAspin Mini
RNA Isolation Kit (GE Healthcare), following the manufacturer’s instructions. Absorbance ratio at 260/280 ηm was used to assess the RNA purity, a ratio of 1.8-2.0
indicating adequate purity.
2.4.2. Labeling and puriﬁcation of targets
RNA samples were prepared and processed according to protocols supplied by the
manufacturer (GE Healthcare). Brieﬂy, cDNAs were synthesized from puriﬁed RNA(1µg) and control bacterial mRNAs. Samples were puriﬁed using the QIAquick
Spin Kit (Qiagen) and concentrated by SpeedVac. Concentrated pellets were used
in a biotinylated-UTP based cRNA synthesis using the CodeLinkTM Expression
Assay Reagent Kit (GE Healthcare). Labeled cRNAs were puriﬁed using RNeasy
Kit (Qiagen) and fragmented with supplied solution at 94o C for 20 min.
2.4.3. Hybridization and washing of the DNA arrays
Fragmented biotin-labeled cRNAs (10µg) were incubated with CodeLinkTM bioarrays under agitation (300 rpm) for 20h. The bioarrays were then washed and incubated with Cy5-streptavidin (30 min). Scanning of the bioarrays was performed
using a GenePix 4000 B Array Scanner (Axon Instruments) and the data were collected using the CodeLinkTM System Software (GE Healthcare), which provides theSeptember 4, 2009

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

Quality Control in Biological Experiments

7

raw data and invalidated data from irregular spots.
3. Results and Discussions
The proposed error’s estimator was applied to two sets of DNA microarrays, both
set up in technical duplicates, in order to illustrate the application.
In Figure 2, the ﬁrst set of microarrays experiments in duplicates is presented,
constituting an example of an acceptable duplicate of microarrays, since the proportion of rejected spots, i.e., number of rejected spots/total number of spots was∼11%. The SVR ﬁtted curve is in blue and the rejected spots (poor quality spots) inred. Notice in Figure 2A that some bias occurs in microarrays experiments. Moreover, in Figure 2B, it is possible to verify that the variance is not constant, i.e., amodest heteroscedasticity is found in the data.

Fig. 2. An example of acceptable reproducibility in duplicated microarrays. A) Fluorescenceintensities of Array 1 versus Array 2 plot of the ﬁrst set of duplicates in log scale; B) Array 1
versus residues of the ﬁrst set of duplicates with Array 1 in log scale. In red, the spots with
T > 10%; in blue, the ﬁtted curve.

Similar characteristics may also be observed in the second set of microarrays(Figure 3). Considering an α = 10%, the proportion of rejected spots was ∼38%,
therefore, this second set has a lower quality than the ﬁrst one. It is interesting thatthe variance is clearly high in the low signal spots (Figure 3B).
Notice that, in both cases, most of the spots marked in red display low expression, which is to be expected, since it is known that low signal spots display higher

September 4, 2009

8

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

A. Fujita et al.

Fig. 3. An example of poor reproducibility in duplicated microarrays. A) Fluorescence intensities
of Array 1 versus Array 2 plot of the second set of duplicates in log scale; B) Array 1 versus
residues of the second set of duplicates, with Array 1 in log scale. In red, the spots with T > 10%;in blue, the ﬁtted curve.

variance than high signal ones [26]. Therefore, it is not possible to distinguish thetrue signal from the noise background signal when the spot signal is low.
By calculating the Pearson’s correlation for both experiments, we have obtained
ρ1 = 0.90 (p-value < 0.01) and ρ2 = 0.93 (p-value < 0.01), respectively. Therefore,
using the Pearson’s correlation results, one may conclude that the replicates are
acceptable, however, the second set of replicates was previously described by ourestimator of error of measure as being unsatisfactory. This illustrates the fact thathigh association (correlation) between duplicates is not the same as reproducibility.
In order to quantify the quality of the spots, several microarray analysis softwares, such as CodeLinkTM System Software [18], ArrayVision v.8.0 (Imaging Research Inc, Ontario, Canada) and TM4 [19] were developed. In most of them theanalysis is based on the ratio between the background signal and the spot signal,and on the spot’s shape (if its shape is well deﬁned) to distinguish acceptable ones.
It is a very important step, which should be performed, but, unfortunately, withthese approaches it is not possible to identify spots with high or low error in measure
derived from technical problems.
SVR does not require a high performance computer to be calculated, i.e., it may
be computed in a personal computer. If millions of measurements per chip become
available, by comparing a pair of technical replicates, we estimate that SVR maySeptember 4, 2009

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

Quality Control in Biological Experiments

9take of the order of a few minutes to compute using a personal computer, i.e., it istotally feasible in practice.
Identiﬁcation of high noise spots using replicates may be an additional criterion
to discard spots which may inﬂuence the following steps in the microarray analysis.
Therefore, sometimes, it is not necessary to discard the whole microarray due toa few genes which present poor measures. One may consider to discard the whole
microarray if the number of rejected spots is higher than a certain ratio (numberof rejected spots / total number of spots), for example, 20%. Otherwise, one may
discard only the rejected spots.
Here, we have illustrated our estimator of error of measure for duplicates, but
for more than duplicates, it may be obtained, in a straightforward manner, bycalculating all di combinations among replicates.
In summary, we propose an interpretable and useful method in order to distinguish acceptable replicates from poor replicates. In addition, we have presented asolution to overcome some well known problems within Dahlberg’s error and modeled the heteroscedasticity present in microarrays. Our illustrative examples arefocused in gene expression data. However, the proposed general method may be
applied to any quantiﬁcation procedure, such as, protein’s quantiﬁcation or RealTime RT-PCR experiments.

3.1. Acknowledgements
This work was supported by RIKEN, Japan; and by Brazilian research agencies(FAPESP, CNPq, FINEP and CAPES).

References[1] Bakay, M., Chen, Y.W., Borup, R., Zhao, P., Nagaraju, K. and Hoﬀman, E.P., Sourcesof variability and eﬀect of experimental approach on expression proﬁling data interpretation, BMC Bioinformatics, 3:4, 2002.
[2] Balagurunathan, Y., Dougherty, E.R., Yidong, C., Bittner, M.L. and Trent, J.M.,
Simulation of cDNA microarrays via a parameterized random signal model, Journalof Biomedical Optics, 7:507-523, 2002.
[3] Bammler, T. et al., Standardizing global gene expression analysis between laboratories across platforms, Nat Methods, 2:351-356, 2005.
[4] Chirgwin, J.M., Prybyla, A.E., Macdonald, R.J. and Rutter, W.J., Isolation of biologically active ribonucleic acid from sources enriched in ribonucleases, Biochem,
18:5294-5299, 1979.
[5] Dahlberg, G., Statistical methods for medical and biological students, Interscience
Publications, New York, 1940.
[6] Dellavia, C., Sforza, C., Orlando, F., Ottolina, P., Pregliasco, F. and Ferrario, V.F.,
Three-dimensional hard tissue palatal size and shape in Down syndrome subjects,
The European Journal of Orthodontics, 29:417-422, 2007.
[7] Draghici, S., Khatri, P., Eklund, A.C. and Szallasi, Z., Reliability and reproducibility
issues in DNA microarray measurements, TRENDS in Genetics, 22:101-109, 2005.
[8] Fujita, A., Sato, J.R., Rodrigues, L.O., Ferreira, C.E. and Sogayar, M.C., Evaluat-

September 4, 2009

10

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

A. Fujita et al.

ing diﬀerent methods of microarray data normalization, BMC Bioinformatics, 7:469,
2006.
[9] Kiryu, H., Oshima, T. and Asai, K., Extracting relations between promoter sequences
and their strengths from microarray data, Bioinformatics, 21:1062-1068, 2004.
[10] Houston, W.J.B., The analysis of errors in orthodontic measurements, American
Journal of Orthodontics, 83:382-390, 1983a.
[11] Houston, W.J.B., Walther’s orthodontic notes, 4th ed. Wright PGS, Bristol., 1983b.
[12] Jarvinen, A.K., Hautaniemi, S., Edgren, H., Auvinen, P., Saarela, J., Kallioniemi,
O.P. and Monni, O., Are data from diﬀerent gene expression microarray platforms
comparable?, Genomics, 83:1164-1168, 2004.
[13] Jenssen, T.K., Langaas, M., Kuo, W.P., Smith-Sorensen, B., Myklebost, O. andHovig, E., Analysis of repeatability in spoted cDNA microarrays, Nucleic Acids Res.,
30:3235-3244, 2002.
[14] Kamoen, A., Dermaut, L. and Verbeek, R., The clinical signiﬁcance of error measurement in the interpretation of treatment results, European Journal of Orthodontics,
23:569-578, 2001.
[15] Kim, B.S., Rha, S.Y., Cho, G.B., and Chung, H.C., Spearman’s footrule as a measureof cDNA microarray reproducibility, Genomics, 84:441-448, 2004.
[16] Naderi, A., Ahmed, A.A., Barbosa-Morais, N.L., Aparicio, S., Brenton, J.D. and Caldas, C., Expression microarray reproducibility is improved by optimising puriﬁcation
steps in RNA ampliﬁcation and labelling, BMC Genomics, 5:9, 2004.
[17] Qiu, J., Sheﬄer, W. and Noble, W.S. Ranking predicted protein structures withsupport vector regression, Proteins, 71:1175-1182, 2008.
[18] Ramakrishnan, R., Dorris, D., Lublinsky, A., Nguyen, A., Domanus, M., Prokhorova,A., Gieser, L., Touma, E., Lockner, R., Tata, M., Zhu, X., Patterson, M., Shippy, R.,
Sendera, T.J. and Mazumder, A., An assessment of Motorola CodeLinkTM microarray performance for gene expression proﬁling applications, Nucleic Acids Research,
30:e30, 2002.
[19] Saeed, A.I., Sharov, V., White, J., Li, J., Liang, W., Bhagabati, N., Braisted, J.,
Klapa, M., Currier, T., Thiagarajan, M., Sturn, A., Snuﬃn, M., Rezantsev, A., Popov,
D., Ryltsov, A., Kostukovich, E., Borisovsky, I., Liu, Z., Vinsavich, A., Trush, V. andQuackenbush, J., TM4: a free, open-source system for microarray data management
and analysis, BioTechniques, 34:374-378, 2003.
[20] Sforza, C., Peretta, R., Grandi, G., Ferronato, G. and Ferrario, V., Three-dimensional
facial morphometry in skeletal Class III patients: A non-invasive study of soft-tissuechanges before and after orthognathic surgery, British Journal of Oral and Maxillofacial Surgery, 45:138-144, 2007.
[21] Smola, A.J. and Schölkopf, B. A tutorial on support vector regression, Statistics andComputing, 14:199-222, 2004.
[22] Sonnesen, L. and Bakke, M., Molar bite force in relation to occlusion, craniofacialdimensions, and head posture in pre-orthodontic children, European Journal of Orthodontics, 27:58-63, 2005.
[23] Vapnik, V. and Lerner, A., A pattern recognition using generalized portrait method,
Automatic and Remote Control, 24:774-780, 1963.
[24] Vapnik, V. and Chervonenkis, A., A note on one class of perceptrons, Automation
and Remote Control, 25, 1964.
[25] Vapnik, V., Statistical learning theory, New York: Wiley, 1998.
[26] Wu, T.D., Large-scale analysis of gene expression proﬁles, Brieﬁngs in Bioinformatics,
3:7-17, 2002.
[27] Yang, Y.H., Buckley, M.J., Dudoit, S. and Speed, T.P., Comparison of methods forSeptember 4, 2009

16:21

WSPC - Proceedings Trim Size: 9.75in x 6.5in

20-Fujita

Quality Control in Biological Experiments

11

image analysis on cDNA microarray data, Journal of Computational and Graphical
Statistics, 11:108-136, 2002.
[28] Yuk, F.L. and Cavalieri, D., Fundamentals of cDNA microarray data analysis, Trends
in Genetics, 19:649-659, 2003.