cDNA arrays
The microarrays used in this paper were produced in house using a Micro Grid II robotic printer (Bio Robotics, Cambridge, UK). These 13 k human cDNA arrays were printed on amino silane coated slides (CMT GAPS, Corning Life Sciences, Corning, NY). For details on the arrays, we refer to: http://www.med.uio.no/dnr/microarray/index.html.
RNA purification
RNA from two cell lines, MT-1 (human breast cell line, kindly provided by Dr. I Fitchner, (Naundorf et al. [23]) and U2OS (human osteosarcoma line, ATCC, Rockville, MD), were used in all experiments. The cell lines were maintained in RPMI media (Bio Whittaker Europe) supplemented with 10% fetal calf serum (PAA Laboratories, Linz, Austria). The cell cultures were grown without antibiotics and tested negative for mycoplasma infection. Total RNA was isolated using GenElute Mammalian Total RNA kit (Sigma, St. Louis, MO) according to the manufacturer's protocol. RNA concentrations were determined by OD260 reading in 50 mM NaOH (GeneQuant, Amersham Pharmacia Biotech AB, Uppsala, Sweden). Total RNA from the respective cell lines were aliquoted to 30 μg total RNA per sample. In addition, an aliquot of RNA from each cell line was diluted to 1 μg/ul total RNA concentration and further down to 0.1 μg/ul.
RNA amplification
RNA amplification was carried out as published by Baugh et al. [5] though slightly modified. The modifications mainly involved a different choice of enzyme in respect to the original protocol. A detailed version of the protocol can be found at the following website: http://www.med.uio.no/dnr/microarray/Mat_and_met_amp.pdf. Briefly, the initial amount of total RNA used for all amplifications was 0.2 μg and all samples went through two rounds of amplification. In a 10 μl reverse transcription reaction, 0.2 μg total RNA was transcribed by priming with a dT/T7-primer. Second strand synthesis was initiated by RNase digestion. The purified double stranded cDNA served as the template for the first round of in vitro aRNA transcription. For the subsequent round of amplification, the first strand cDNA was primed with random hexamers. The synthesis of the second strand was initiated by annealing a dT/T7-primer to the aRNA:cDNA heteroduplex with partially digested aRNA. After purification, a second round of in vitro transcription followed.
RNA labeling
Three μg of aRNA and 30 μg of total RNA was labeled with Cy3- or Cy5-labeled dCTP (Amersham Pharmacia Biotech AB) during reverse transcription. The reaction mixture consisted of random hexamers (16 μg) or anchored oligo-dT 20-mer primers (4 μg) depending on whether aRNA or total RNA was labeled. Further, 40 U RNAsin (Promega, Madison, WI), 1st strand buffer, 0.1 M DTT, 0.5 mM of dATP, dCTP, dGTP and 0.2 mM dTTP was added and the mix was incubated for 5 min at 65°C in a water bath. The tube was then transferred to a heating block (42°C), and 4 μl of either fluorophore was added to the respective tubes followed by 400 U SuperScript II (Invitrogen, Groningen, The Netherlands). After 60 min the reaction was stopped by adding 5 μl 0.5 M EDTA (pH = 8.0). Residual RNA was hydrolyzed with 10 μl 1 M NaOH and incubated at 65°C for either 60 or 15 min depending on whether total RNA or aRNA was reverse transcribed. Twenty-five μl 1 M Tris-HCl (pH 7.5) was added to neutralize the mixture. Micro Bio-Spin columns, P-6 (Bio-Rad Laboratories, Hercules, CA) were used to remove unincorporated dye. Labeled Cy3- and Cy5-cDNA was respectively diluted with TE-buffer (pH 7.5) before Microcon YM-30 columns (Amicon, Millipore Corporation, Bedford, MA) were applied to concentrate the samples to 17–20 μl.
Hybridization and scanning
The hybridization volume of 40 μl consisted of 17 μl of each of the labeled probes, 3.5 × SSC (pH = 7.5), 0.3% SDS, 1.25 × Denhardt's solution, 4 μg yeast tRNA, and 4 μg BSA. Additionally, 16 μg poly A (Amersham Pharmacia Biotech AB) was used if the probe was based on total RNA. The final mix was heated for 2 min at 100°C and spun down for 10 min at 13 K before it was applied on a microarray slide under the LifterSlip (Erie Scientific Company, Portsmouth, NH). The slide was then placed in an ArrayIT hybridization chamber (Telechem, Sunnyvale, CA) and submersed into a water bath for over night hybridization at 65°C. Prior to scanning, the coverslip was removed in a solution of 2 X SSC and 0.1% SDS, and the slide was washed in the following solutions for 5 min at room temperature; 1 X SSC, 0.2 X SSC and 0.05 X SSC. The slide was dried by centrifugation. Scanning was performed with a ScanARRAY 4000 (Packard Biosciences, Biochip Technologies LLC, Meriden, CT) scanner, and data was acquired from the images using GenePix Pro 3.0 software (Axon Instruments Inc., Union City, CA). A background-corrected intensity for each spot was formed by subtracting the mean of the pixels in the local background from the mean of the pixels in the spot.
Description of experiments
To avoid variability due to total RNA preparation, we pooled purified total RNA into a single large sample for each of the two RNA sources, the MT-1 and U2OS cell lines. The RNAs used in the experiments were based on these two initial samples. Aliquots from the initial sample represented a batch and there were six batches per cell line. Two of these batches were directly labeled in duplicated experiments using a dye swap strategy, while the remaining 4 diluted batches were used as starting material for independent, parallel amplifications. Corresponding pairs of aRNA, aRNA (MT-1) and aRNA (U2OS), were labeled and co-hybridized in duplicated experiments using (again) a dye swap strategy. Thus, in total there were 4 data sets based on non-amplified RNA and 8 based on amplified RNA. The microarray slides were all from the same print batch. Since it is impossible to hybridize non-amplified and amplified material on the same array, we chose to use two different cell lines on each microarray.
Data preparation
From the data for each microarray we first manually removed technically flawed spots. Then spots automatically flagged by the GenePix software were removed. Finally, all spots where the spot intensity (uncorrected foreground intensity) was lower than twice the background intensity in any of the two channels were removed (i.e. ((spot-background)/background) < 1). An overview of the number of filtered genes is given in Tables 1 and 2.
Systematic errors were corrected by normalizing the data using the locally weighted scatterplot smoother, lowess, as described in Yang et al. [24] and implemented in the sma package in the statistical language R. This routine normalized the ratios for each spot on each array by subtracting a value s, depending on the value of the log2-transformed total intensity of the spot on the array. In addition we normalized log2-transformed intensities by adding the value s/2 to the log2-transformed intensity of the Cy3 (green) channel and subtracting the value s/2 from the log2-transformed intensity of the Cy5 (red) channel. This normalization of intensities was done because the ANOVA model is based on log2-transformed intensities.
Analysis of variance modeling
The first step in the ANOVA modeling is to evaluate the sources of variability present, and their interactions. The next step is to construct the model and experimental design accordingly. We identified the following 6 main effects or factors for the design:
- We used 12 different microarrays (A).
- There were 2 cell lines (C).
- Each cell line was divided into 6 batches, giving a total of 12 batches (B).
- The material was labeled red (Cy5) or green (Cy3), so the dye (D) effect had 2 levels.
- The material was amplified on 8 of the microarrays and not amplified on the other 4, i.e. two levels for the protocol (P) effect.
- A total of 10759 genes (G) remained after filtering (these genes were present on at least one microarray).
Figure 1 presents an overview of the experimental design. For a particular cell line, the material was divided into 6 batches. These batches were split in two and labeled Cy5 (red) or Cy3 (green). The design is not balanced, since we chose to use twice as many microarrays with amplified material, but we wanted to be able to handle the possibility that each amplification event induced a larger variability between the arrays with amplified material than between the arrays with the non-amplified material.
Statistical model for the ANOVA analyses
To investigate the different sources of variability, we set up an ANOVA-based statistical model. Related models are found in Kerr et al. [25], Wolfinger et al. [26] and Jin et al. [27]. Let log2 y
ijklmn
denote the log2-transform of the normalized measured intensity on array i, for cell line j, for batch k, using dye l and protocol m for gene n. We explain the log2-transformed intensity by the following model:
log2 y
ijklmn
= μ + A
i
+ C
j
+ B
k
+ D
l
+ P
m
+ G
n
+AD
il
+ CP
jm
+ GP
nm
+ CG
jn
+ AG
in
+ DG
ln
+ CGP
jnm
+ ε
ijklmn
,
with i = 1,...,12, j = 1,2, k = 1,...,12, l = 1,2, m = 1,2, n = 1,...,10759, where μ is the overall mean, A
i
is the overall array effect on array i, C
j
is the overall cell line effect of cell-line j, B
k
is the overall batch effect of batch k, D
l
is the overall dye effect of dye l, P
m
is the overall effect of the amplification protocol and G
n
is the overall gene effect of gene n. Furthermore, AD
il
is the interaction between array and dye, so if this effect is significant, the dye has different effect for different arrays. CP
jm
models different effect of amplification on different cell lines. GP
nm
is of significant size if the amplification protocol has a different effect on different genes and CG
jn
represents the gene-specific cell line effect. This is usually an effect of interest e.g. if we want to measure different expression for the two cell lines for each gene. AG
in
is the spot effect, it models variation in the amount of cDNA deposited in each spot during the printing process. DG
ln
represents the gene-specific dye effect. CGP
jmn
is the cell line specific GP
nm
effect, and is our interaction of primary interest.
Only the cell line effect (C
j
), the dye effect (D
l
), the protocol effect (P
m
) and the cell line protocol interaction (CP
jm
) were modeled to be fixed effects (the first three have only two levels and the last is a combination of two fixed effects). CP
jm
is confounded with C
j
and P
m
, so we only estimated CP22. All other effects were modeled as independent normally distributed random effects with zero mean and constant standard deviation. This means that A
i
~ N(0,
), B
k
~ N(0,
), and similarly for all the other random effects. ε
ijklmn
is the modeling and measurement error, and it was also assumed to be normally distributed, ε
ijklmn
~ N(0,
). This is a mixed-effects model, since some of the effects are fixed and others are random.
We were interested in the preservation of the ratios after amplification and used the model for the log2-transformed intensity, log2 y
ijklmn
, to set up an expression for the ratio of cell line 1 and cell line 2 for each spot on each array. Considering a randomly chosen spot, we wrote the observation of cell line 1 (j = 1) as log2 y
ijklmn
and cell line 2 (j = 2) as log2 yijk'l'mn(the same array, i, was used, the same protocol, m, was used, it was the same gene, n, but batch, k and k', and dye, l and l', were different). The effects that were equal for each spot cancelled when looking at the ratio:
We wanted to see how the overall signal was affected by an amplification protocol that possibly did not fully preserve the ratios. The signal we are usually interested in is the gene expression using a specific cell type for a specific gene, i.e. CG
jn
. The Signal to Noise Ratio (SNR) is a valuable measure for answering this question. The SNR is defined as signal divided by noise. CG
jn
is the signal and the other effects are noise. Using the expression for the ratios, we see that the variance of the ratio of interest is
(the variance for the fixed effects are VAR(C
j
) = VAR(D
l
) = VAR(CP
jm
) = 0) so the Signal to Noise Ratio is
By comparing the Signal to Noise Ratio with and without the CGP
jmn
interaction present in the model, we were able to see how large impact amplification had on the signal. Hence, we estimated the parameters in the ANOVA model, and used the parameters estimated for the random effects to form the SNR to be investigated.
Due to the large dimensions in the mixed-effects model, we estimated the parameters of interest (the fixed effects and the variances of the random effects) using Gibbs-sampling. Gibbs-sampling is an iterative simulation method, which was run on the computer to provide our parameter estimates. We refer to the web page http://www.nr.no/documents/samba/research_areas/SMBI/NAR/ for details.
Hypothesis testing
Preservation of ratios can be addressed by performing hypothesis testing by comparing ratios between the two cell lines for the amplified material with ratios between the two cell lines for the non-amplified material for each gene. We have used this strategy as a complement to our preferred ANOVA-based procedure. Let ν
n
(1) be the expected log2-ratio between the two cell lines for the non-amplified material for gene n, and let ν
n
(2) be the expected log2-ratio between the two cell-lines for the amplified material for gene n. We wanted to test the following hypothesis for each gene n.
H0n: ν
n
(1) = ν
n
(2) vs H1n: ν
n
(1) ≠ ν
n
(2)
Many methods for performing hypothesis testing exist, and we chose to use a two-sample t-test for each gene separately. Thus, we estimated the variability for each gene in each group of samples separately, amplified and non-amplified (allowing the variability for the amplified samples to be different from the variability for the non-amplified sample). The Smith-Welch-Satterthwaite approximation to the degrees of freedom in the two-sample t-test was used, see any introductory statistical text e.g. Walpole, Myers and Myers [28], chapter 10.8.