Outliers involving the Poly(A) effect among highly-expressed genes in microarrays

Hybridization to the microarray

Fig. 1 shows background-subtracted images of a microarray that has been hybridized with first strand cDNA derived from skeletal muscle mRNA. Also shown in Fig. 1 are the background-subtracted images of the array after hybridization with end-labeled oligo(dT), performed in order to measure accessible PolyA in the spots.

Search for artifacts from the PolyA effect

As shown in the scatterplot provided in Fig. 2, there is only a weak correlation between the hybridization signals of spots with first strand cDNA from rat muscle, versus hybridization of the spots with oligo(dT). Among the 567 microarray spots having muscle gene-expression signals in the range that could be detected (as defined in the Methods section), the Pearson correlation between the muscle signals and the signals obtained using the oligo(dT) probe is only 0.46, and the correlation is even weaker (0.37) when the correlation calculation is restricted to microarray spots having a muscle signal greater than the median value of 7682 phosphorimager units. We therefore conclude that a polyA effect is not readily apparent in these data, presumably because we had added unlabeled polyA to the hybridization mixture as a blocking agent.

Spots at random locations throughout the array did not bind significantly to the oligo(dT) probe, almost all of which correspond to spots at which there was also no binding of muscle cDNA. Such spots do not appear to represent defects in the manufacturing of the array filter at which no DNA was deposited, because at the end of the experiment, we looked for DNA in those spots directly by staining them with the fluorescent nucleic acid dye acridine orange, followed by imaging of the spots with a fluorimager. Spots that gave no signal with the oligo(dT) probe nevertheless had a corresponding fluorescence signal, except for the deliberately vacant positions in the right hand side of each of the microarray matrices shown in Fig. 1. Furthermore, there was one spot at which binding to the oligo(dT) probe was undetectable, but at which there was nevertheless a very strong binding to the muscle probe. That spot (field 1d, column 11, row 21 in Fig. 1 and position 14511,0 in Fig. 2) corresponds to the enteric smooth muscle form of gamma actin (Accession number T60048), which according to GenBank's annotation for that accession number is a possible reversed clone because polyT (or polyA in its complementary strand) was not found within it.

Observation of outlier microarray spots

Although there is generally little correlation between the hybridization values obtained with oligo(dT) and muscle cDNA, this may not be true for spots that show the greatest binding to oligo(dT). For the spots showing oligo(dT) hybridization values greater than 12000 (phosphorimager units), the corresponding value of the muscle cDNA hybridization is always 24000 or more (phosphorimager units). This result was obtained for the same eight microarray spots, when microarray hybridizations were performed using muscle mRNA from both fa/fa and Fa/Fa Zucker rats. In contrast, spots with oligo(dT) hybridization values less than the threshold value of 12000 exhibit corresponding values for the muscle probe hybridization that may be anywhere between zero and the upper range of values.

The finding is not likely to be a chance happenstance for the following reason. If the null hypothesis is that the muscle-probe hybridization values for spots with polyA values > 12000 have the same statistical distribution as spots with polyA values < 12000, then the probability that all eight of them have values greater than the median is (1/2)⁸ = 0.0039, and the probability that all eight of them have a value greater than 24000 is less than (52/567)⁸ which is less than 10^-8, because only 52 of the 567 detectable spots have values greater than 24000. Thus, we reject the null hypothesis and conclude that the statistical distribution of muscle-probe values among spots having polyA values greater than the threshold of 12000 is significantly different than that for the remainder of the spots.

The question then arises as to whether there is anything unusual about the clones used to generate the DNA for those eight outlier microarray spots that a priori would have allowed us to flag them as potential concerns? As now described, distinguishing features of the clones might be that (1) they have unusually short inserts or (2) they may contain an unusually large fraction of bases that might be expected to bind to oligo(dT), namely, the bases that lie within a sequence run of AAAAAA... .

Relation between clone sequence and oligo(dT) probe hybridization signal

A search of GenBank reveals that only partial sequences are available for the clones corresponding to the eight outlier spots. One of the clones has a partial sequence containing a polyA stretch consisting of 41 bases (Accession H94857 (3') and H94912 (5') = GCN5 general control of amino-acid synthesis 5-like 1 = spot 1e6,28 in Fig. 1 = point 33211,15102 in Fig. 2). Five of the others have stretches of A or T between 15 and 26 bases, which are not unusual (Accession R52548 (3') and R52604 (5') = Superoxide dismutase 1 = spot 1c1,21 in Fig. 1 = point 24216,12384 in Fig. 2; AA460830 (3') and AA461132 (5') = RNA polymerase II polypeptide J = spot 1b11,3 in Fig. 1 = point 48345,12055 in Fig. 2; AA434115 (3') and AA434048 (5') = Cartilage glycoprotein 39 = spot 1b7,20 in Fig. 1 = point 30221,14554 in Fig. 2; AA490256 (3') and AA490356 (5') = Guanine nucleotide-binding protein G(K) alpha subunit = spot 1f,7,19 in Fig. 1 = point 24195,14273 in Fig. 2; and H93118 (3') and H93246 (5') = EST = spot 1c12,26 in Fig. 1 = point 42177,15785 in Fig. 2). The other two clones contain stretches of A no longer than 6 bases (Accession AA456352 (5') and AA454703 (3') = DEAD/H Asp-Glu-Ala-Asp/His box peptide 38 = spot 1f4,14 in Fig. 1 and point 32545,14353 in Fig. 2; and N92646 (3') and N99582 (5') = granulocyte-macrophage colony stimulating factor receptor 1 = spot 1f5,21 in Fig. 1 = point 41079,15760 in Fig. 2). However, N92646 contains an A-rich stretch of 30 base pairs with 23 A bases, and AA454703 contains more than one stretch of 6-A bases. Because stable hybridization duplexes have been reported with oligonucleotides as short as 7 bases [5], and because such stable duplexes may be formed even if there are base-pair mismatches [6], it is conceivable that the oligo(dT) probe could have hybridized to the A-rich regions in the 3' or 5' ends of the DNA in each the outlier spots.

The oligo(dT) probe could also have hybridized to the unknown sequences between the 3' and 5' ends of these clone inserts, so we sought their complete sequences. For purposes of comparison, we also sought the complete sequences of all clone inserts for which the muscle cDNA gave a hybridization signal in the range exhibited by the outlier clones (>24000 phosphorimager units). We sought the full sequences by using the clustering tool "IMAGEne" available at the IMAGE Consortium web page http://image.llnl.gov. This tool attempts to match the 5' and 3' partial sequences of a clone insert with the sequences of all other IMAGE clone inserts, as well as with all RefSeq reference sequences [7]. If matches are found, the entire clone insert sequence may then be constructed by piecing together overlapping matches of the contiguous sequences. Among the 52 clones for which the muscle cDNA gave a hybridization signal in the range exhibited by the outlier clones (>24000 phosphorimager units), it was possible to construct the entire insert sequence of 18 clones (see Additional File 1). Three of these clones are among the outlier group for which the polyA signal is greater than 12000 phosphorimager units.

We then noted two features of the sequence data that may explain in part the large signals that the outlier microarray spots exhibited when hybridized to the oligo(dT) probe. The first feature is that the hybridization signal is generally inversely proportional to sequence length, and two of the outlier clones have unusually short lengths. As seen in Fig. 3, the two outlier clones having insert lengths of less than 400 base pairs give the strongest hybridization signals. Because each microarray spot has the same total amount of DNA (5 ng), the microarray spots associated with these clones must have an unusually large number of short, densely deposited target molecules, and should therefore produce a proportionately large signal. The third outlier spot had an oligo(dT) hybridization signal only barely larger than the threshold of 12000 phosphorimager units and has a somewhat longer insert length (622 base pairs). The second noteworthy feature of the sequence data is that outlier spots may contain an unusually large fraction of bases that might be expected to bind to oligo(dT), namely, the bases that lie within a sequence run of seven or more As (AAAAAAA). As seen in Fig. 4, for two of the outlier clones, the fraction of such bases have relatively high values of 0.028 and 0.040. The third outlier spot has an oligo(dT) hybridization signal only barely larger than threshold of 12000 phosphorimager units and has a more typical fraction of 0.020.

Housing and care of animals

Female lean (Fa/Fa) and obese (fa/fa) Zucker rats, 5 wk of age, were purchased from the Animal Model CORE Facility of the University of California at Davis. They were housed in the Animal Resource Center of the University of Texas at Austin under standard laboratory conditions. All procedures were approved by the University of Texas at Austin Animal Use Committee.

Tissue samples

Rats were quickly anesthetized by intravenous injection of sodium pentobarbital (50 mg/kg). Red quadriceps muscles from one leg of each animal were removed and freeze clamped as rapidly as possible with tongs that had been cooled in liquid N₂. Samples were stored at -70°C for later use.

Extraction of total and messenger RNA

Tissue samples were homogenized in an acid-phenol reagent, with the volume of the tissue not exceeding 10 percent of the volume of the reagent used. Total RNA was then obtained from the homogenate by the procedure recommended by the reagent's supplier (TRI-reagent, Sigma #T9424). This was followed by selection of poly(A) mRNA from the total RNA by hybridizing the RNA with oligo(dT) magnetic microparticles, then isolating the mRNA magnetically (mRNA Isolation Kit, Miltenyi Biotec #751-02). The spectrophotometric 260/280 ratios for mRNA from the obese and lean tissues were 1.95 and 1.99, respectively, which were sufficiently close to 2.0 that an additional round of oligo(dT) selection was not performed.

Preparation of first strand cDNA

Complementary DNA was made using 200 ng of purified, oligo(dT) primed mRNA and SuperScript II RNase H- reverse transcriptase (BRL Life Technologies #18064-014). The cDNA was made radioactive by incorporation of ³³P-dCTP (ICN #58201) using the procedure described at http://www.resgen.com/gf200pro.html, except that 33 mM of dATP, dTTP, and dGTP (Pharmacia #27-2035-01) and 1 U/μl of a ribonuclease inhibitor (Ambion #2690) were added. Unincorporated nucleotides were separated from the labeled cDNA using a push column (Stratagene #400701).

Hybridization to array membrane

A high density DNA array membrane was used to measure the level of expression of each of 5,184 genes (Research Genetics, GENEFILTER #GF200, Lot #980611-21). Hybridization of the membrane with radioactive cDNA was performed by the procedures described at http://www.resgen.com/gf200pro.html, except that the last room temperature wash was performed in 2X SSC. The hybridization solution included 0.5 μg/ml poly(dA) (Research Genetics #POLYA.GF) and 0.5 μg/ml cot1 DNA (BRL #15279-011) as blocking agents. Stripping the filter, in order to reprobe it, was done by placing the filter in a boiling 1% SDS solution, as recommended by the manufacturer. The stripped filter was used to expose a phosphorimager plate, which was then scanned. The resulting image had < 0.005 remaining from the image that was obtained before stripping. We obtained nearly identical microarray hybridizations using cDNA prepared from skeletal muscle mRNA from lean (Fa/Fa) versus obese (fa/fa) Zucker rats at an age of 6 weeks (the age at which phenotypic differences between fa/fa and Fa/Fa are thought to begin). Only one gene was possibly expressed differentially (the sarcomeric isoform of the mitochondrial creatine kinase gene, sMtCK, which appeared to have higher expression in the lean muscle). We confirmed the sMtCK result with a Northern blot, but did not do so for any of the other genes, which may be considered a limitation of our experiment.

Quantifying the hybridization results

Radioactivity in hybridized filters was imaged using a phosphorimager (Model 445 SI, Molecular Dynamics). The resulting image files were then processed using custom software, as shown in Fig. 5. Each of the membrane's sixteen 12 × 30-spot arrays was extracted from the original image file, along with a border of pixels corresponding to one spot width. A background value was then estimated for each pixel in each of these images, as follows. If a pixel was situated in the region between or adjacent to DNA spots, and if the spatial derivative of the image at that pixel indicated that it was not part of the overlap region between two intense spots, that pixel was defined to be a background pixel. If a pixel was situated within the region corresponding to a spot of DNA, or if it was situated in the region between two overlapping intense spots, the background value for that pixel was set equal to the value of the nearest background pixel, as defined above. A background-subtracted image of the radioactivity was then obtained by subtracting the background value for each pixel from the original image.

The intensity of spots in an image is a function of how long the phosphorimager plate was exposed, as well as factors such as the efficiency of probe labeling. Therefore, in order to be able to compare background-subtracted images of the same 12 × 30 spot array for successive hybridizations, normalization was performed as follows. The image for one of the hybridizations was selected to be the reference, and corresponding pixels in the other image were normalized by substituting them into a quadratic polynomial normalization function (of the pixel value), the coefficients of which were estimated by singular value decomposition (SVD) [10]. This quadratic model allows for some compensation in the event that the response of the phosphorimager changes between measurements. We found that the fitted second order coefficients were always very close to zero, and the constant offset coefficients were also very close to zero. Thus, normalization of background-subtracted images of the same array for different hybridizations was accomplished essentially by multiplying all pixels in the non-reference image by a constant scale-factor, the value of which had been estimated by SVD. Note that the phosphorimager data have units related to the voltage of its phototube, which are not related in any obvious way to nucleic acid concentrations or amounts. Some investigators rescale the data such that the average signal per microarray spot equals 1, but we did not do so.

The distance between spots on the array is only 0.75 mm, and as a consequence, many of the adjacent intense spots overlap one another. When integrated intensity within a specified region about the center of the spot is used to represent the magnitude of hybridization, it is then useful to model the image as the sum of superimposed two-dimensional spot functions. This was done using two-dimensional gaussian functions to represent the spots, as shown in Fig. 5, with the location, amplitude, and standard deviation coefficients all estimated automatically from the background-subtracted image in two dimensions, using the Levenberg-Marquardt method [10]. The similarity between background-subtracted and fitted spot images (Fig. 5b vs. 5c) indicates that much of the spread of the intense spots appears to be gaussian, which might be attributed to the gaussian shape of the laser beam that scans the phosphorplate. When the spot intensity was close to the noisy background level, this model cannot be used to fit the data, due to near singularity of the matrix equations that had to be solved to perform the fitting. In that case, the method that we used to estimate the intensity of each spot was to sum the nine adjacent pixels in the center of each spot of the background-subtracted image. The correlation between the values so obtained and the corresponding value of the integrated spot intensity from successful gaussian curve fitting was 0.99.

A histogram was constructed from the logarithm of the values of all the spot values. As observed in [11], the histogram for our muscle data was bimodal, consisting of a gaussian-like distribution of low-intensity spots, overlapping an adjacent distribution of moderate and high-intensity spots. The transition between these two distributions occurred within a clearly recognizable range of values, 1550 to 1800 phosphorimager units. We therefore followed the conservative practice recommended in [11] by considering all spot values less than a value of 1800 to define undetectable hybridizations. According to this criterion, 567 of the spots were detectable.

Search for promiscuous hybridization to polyA segments

To investigate the polyA effect as a potential artifact, we performed the following experiment, using essentially the approach described in [4]. Oligo(dT) (10–20 mer mixture, Research Genetics cat. # POLYT.GF) was end labeled with T4 kinase (Life Technologies #10476-018) and (γ-³³P)ATP (ICN #58000). The array filter was then hybridized with this probe in Hybrisol I (Oncor #S4040) for 10 hours at 42°C, then washed twice at room temperature for two minutes in 6X SSC and 0.1% SDS. It was then used to expose a phosphorplate for 12 hours, which was scanned using a phosphorimager (Model 445 SI, Molecular Dynamics) to identify the array spots having long stretches of poly(dA). Quantifying of the spot intensities was performed as indicated above.

Staining of the microarray filter with a fluorescent nucleic acid indicator

After all hybridizations of the microarray filter had been performed, we stained the filter with a fluorescent dye that is specific for nucleic acids. We did so in order to determine whether the absence of hybridization signals for some microarray spots was due to the absence of any DNA there, which would indicate a defect in the manufacturing of the microarray at those spots. The filter was stained with 1 μg/ml acridine orange (Sigma A6014) for 10 minutes at room temperature, rinsed in distilled water, then imaged with a fluorimager (Molecular Dynamics). The dye was excited by the 488 nm argon laser line of the fluorimager, and emission was detected with a narrowband filter centered at 530 nm.