Identification and correction of previously unreported spatial phenomena using raw Illumina BeadArray data
© Smith et al. 2010
Received: 15 December 2009
Accepted: 27 April 2010
Published: 27 April 2010
Skip to main content
© Smith et al. 2010
Received: 15 December 2009
Accepted: 27 April 2010
Published: 27 April 2010
A key stage for all microarray analyses is the extraction of feature-intensities from an image. If this step goes wrong, then subsequent preprocessing and processing stages will stand little chance of rectifying the matter. Illumina employ random construction of their BeadArrays, making feature-intensity extraction even more important for the Illumina platform than for other technologies. In this paper we show that using raw Illumina data it is possible to identify, control, and perhaps correct for a range of spatial-related phenomena that affect feature-intensity extraction.
We note that feature intensities can be unnaturally high when in the proximity of a number of phenomena relating either to the images themselves or to the layout of the beads on an array. Additionally we note that beads neighbour beads of the same type more often than one might expect, which may cause concern in some models of hybridization. We highlight issues in the identification of a bead's location, and in particular how this both affects and is affected by its intensity. Finally we show that beads can be wrongly identified in the image on either a local or array-wide scale, with obvious implications for data quality.
The image processing issues identified will often pass unnoticed by an analysis of the standard data returned from an experiment. We detail some simple diagnostics that can be implemented to identify problems of this nature, and outline approaches to correcting for such problems. These approaches require access to the raw data from the arrays, not just the summarized data usually returned, making the acquisition of such raw data highly desirable.
A key stage for all microarray analyses is the extraction of feature-intensities from an image. If this step goes wrong, then the experiment will certainly be compromised. Thus, much research has gone into the tasks of automatic identification of both the features on arrays and the grids on which they lie [1–4]. Illumina BeadArrays differ from other types of microarray in their construction, and have had special attention paid to them as a consequence [5, 6].
Illumina employ random construction of their arrays on a hexagonal grid . As a consequence, each probe will occur a random number of times on the array and in random locations. This means that the beads have to be found twice with Illumina BeadArrays (thrice with two-colour versions): once by Illumina to identify the type of bead present at a location, and once (or twice) by the user to quantify intensities after hybridization. These steps are then even more important for the Illumina platform than they are for other technologies.
While Illumina's software can report the raw bead-level data from an array (i.e. the location and intensity for every bead), more typically only summarized data are produced. We have previously shown that there are advantages to working with the bead-level data , not least the abilities to calculate covariances in two colour platforms  and to identify and correct spatial artefacts . In this paper we show that using the bead-level data it is possible to identify, control, and perhaps correct for a range of other phenomena related to the locations of the beads. In particular, we will consider a number of potential effects that are suggested by previous studies.
For example, it has been shown  for microarrays that the depletion of target molecules, due to those molecules hybridizing to probes, can affect the physical chemistry to an extent that the registered intensity will also be affected. Such effects can occur locally, and the nature of the Illumina array means beads hybridizing to the same target may be proximal, in which case concerns about local depletion may be raised. A previous study  showed anecdotal evidence that clusters of beads that failed to be decoded were spatially associated with regions of outliers, suggesting that such regions might be beneficially excluded.
Other theoretical aspects of the technology require investigation. The background intensity for arrays is calculated in a non-robust manner that might lead to extreme and biased values. Additionally, the high-density structure of the array may lead to situations where the intensity of a bead may influence that of its neighbours.
There are three key steps in moving from an image of a BeadArray to a set of intensities. First the bead locations must be identified, second the beads must be mapped to a bead-type, and third, an intensity value extracted. Here we investigate each of these steps, showing how each can bias the final intensity, but showing also how such errors can be identified and corrected when working with bead-level data. We illustrate this using two bead-level data sets (one single-colour expression, one two-colour genotyping/DNA copy number).
The copy number data we use consist of two arrays from a previously published data set, full details of which are provided by Curtis et al . The arrays we consider have IDs 4127130020 and 4127130188, to which were hybridized a tumour sample and HapMap CEPH NA10851 cell line respectively. The expression data set consists of two chips (IDs 4343238066 and 4343238080) that have been hybridized with two common commercial reference RNA sources (Stratagene Universal Human Reference RNA, and the Ambion Human Brain Reference Total RNA). Since public repositories are not designed to deal with the raw Illumina data that we use here, these data are available to download from our website .
The raw data we consider to be the .tif image file arising from the scanned array, the .txt file containing the locations of all decoded beads (decoding is discussed in the next section) as well as the associated (background-corrected) intensities, the .locs file that additionally contains the locations of the beads that were not decoded, and the .sdf file that contains details of the chip structure. To reduce download sizes the .txt and .locs file have been compressed using the package BeadDataPackR . Instructions for obtaining bead-level data from the scanner can also be found on our website.
All analyses were conducted in R, primarily through the beadarray package . Where we introduce functionality that is not part of beadarray or another BioConductor package, we make this available through an R script containing the additional functions (Additional File 1), and a vignette that will reproduce the figures and tables in this manuscript (Additional File 2).
The array having been scanned, we can consider Illumina's preprocessing as consisting of three distinct steps: registration of beads, mapping of bead IDs, and feature-intensity extraction.
The registration process is described in full by Galinsky , but can be summarized thus: First, a bead detection algorithm is used to locate bright beads, then based on these positions, the locations of the remaining beads are interpolated. Finally the grid is shifted to ensure it is centred over the array. Due to the random construction, in addition to finding the beads it is necessary to match the beads to their bead-types. Illumina do this, using a hybridization-based procedure (termed decoding) to identify the probe type attached to each bead .
The feature-intensity extraction consists of the calculation of background values, the 'sharpening' of the image via a simple filter that compresses observed fluorescence into tighter regions around the features  and then the calculation of foreground values. Background values are simply subtracted from foreground values to return the final intensity.
where P x and P y are the bead-centre coordinates (in pixels), and G x and G y are the grid locations. The statistic gives a measure of departure from the grid. Note that for two-colour platforms, the difference between red and green locations will also be an indicator.
Further, we consider the possibility that the mapping of bead identities to the beads can go wrong. This will either occur for a subset of beads, in which case we may detect it using our measure of departure from the grid, or it will occur for all beads, in which case we must remap the grid positions for each bead ( etc.) and use the bead identities associated with these new grid positions. The shift Δ x can be determined manually, or automatically by finding the shift that minimizes the mean within bead-type variance for the affected segment and maximizes correlation between replicate segments.
A large variety of algorithms have been employed to extract feature-intensity values from scanned microarray images, and their relative merits have been extensively reviewed . Here we consider only Illumina's foreground algorithm  which begins with the bead-centre positions that were identified during the registration step. The foreground value is calculated as a weighted average of sharpened intensities from a 4 × 4 pixel square located about the bead-centre. The centre four pixels of this square always take a maximum weighting, but the weights of the remaining pixels are determined by the fractional part of the bead-centre coordinates as illustrated in Additional File 3.
This process relies on being able to identify the centre of a bead to a resolution of a fraction of a pixel. Illumina report the bead locations to between 2 and 4 decimal places of a pixel in the .txt file, whilst they are stored as single precision floating point values (7-8 decimal places) within the .locs file.
Background values are calculated from a 17 × 17 square of pixels located about the bead-centre in the non-sharpened image . Thus the background being calculated is local rather than global, but not bead-specific. Within this square, the mean of the five lowest intensities is taken to be the background value for that bead. The mean is not a robust summary of such extreme values, and we consider alternatives through new functionality in the beadarray package.
Note that a 17 × 17 square could contain 12 bead-centres (see Additional File 4) and so it is inevitable that beads will share some or all of their five lowest background pixels (or conversely a pixel may contribute to the background of up to 12 beads). It is also clear that many of the pixels in that 17 × 17 square will be contributing to foreground calculations. We mask the 4 × 4 pixel square around each bead-centre that contributes positively to the foreground calculation in order to calculate the number of "true background" pixels a bead actually has out of the 289 considered.
The separate foreground and background values are not returned by Illumina's software, so we have used the beadarray package to obtain these values.
From the details of the foreground calculation, it is clear that the reported bead-intensities are influenced by the precise bead-location. However the registration of bead locations is partially dependent on bead intensities and will influence the precise location of the bead. While in isolation, each of these steps appears sound, the potential for some 'feedback' in this loop (intensity affects location affects intensity) is of concern.
To investigate the magnitude of the bead location effect, we perform simulations of idealized beads. Each bead on the array we consider to be a sphere evenly covered in probes. This we assume results in a smoothly fluorescing feature when the bead is scanned. Then, we simulate the digitization of the image into discrete pixels by creating a pixel grid, integrating the fluorescence within each square of the grid, and rounding to an integer value. We then systematically move the grid, a fraction of a pixel at a time, and calculate the foreground intensity using the known bead-centre location.
We use the fractional part of the bead-centre locations and the intensities obtained from Illumina's .txt file to investigate the patterns of association in real data. We break a theoretical pixel into 100 × 100 bins and, for each bin, plot the mean of all log-intensities for beads whose fractional bead-centre locations coincide with the bin.
Where we identify an association, we investigate the nature of any causality by a) recalculating intensities with location-independent weights, and b) by averaging not the log-intensities but the residual log-intensities once the average for that bead-type has been removed. Since beads of the same type should be showing similar intensities, they should be similarly affected by the bead registration algorithm, and so this second step tests the possibility of log-intensity affecting location.
Clusters of non-decoded beads are identified by finding the locations of all non-decoded beads from the .locs files and mapping the neighbours of all beads. We assume that any such cluster of beads will contain a non-decoded bead with six non-decoded neighbours, and each such bead is used as a seed to identify maximal networks of neighbouring non-decoded beads (by use of an invasion process). Clusters we define as being any such network of 50 or more beads, and these are then expanded further to identify successfully decoded beads in the vicinity of the cluster.
To identify bright beads that might interfere with their neighbours, we use the EBImage package (version 3.3.1)  to segment the image and then select bright features with a mass of at least 40 pixels and a high degree of circularity. Beads with aberrant background calculations are identified using the readTIFF function in beadarray which allows us to interrogate the pixels surrounding a bead.
Identifying beads that neighbour beads of the same type is trivial using the neighbours matrix. In an idealized array containing 1, 000, 000 beads from 50, 000 bead-types where each bead has six neighbours we would expect to have 60 such pairs. This calculation ignores some of the structure of the array (most obviously that not all beads have six neighbours), and we investigate the numbers that we would expect for our arrays under a truly random construction as follows: We maintain the neighbours matrix for an array, but permute all bead identities of beads that were successfully decoded, tally the number of neighbouring pairs of the same type, and repeat.
To assess the impact of the identified phenomena we use the beadarray package along with the scripts we provide here to identify affected beads and remove them prior to summarization. Analyses follow the scheme illustrated in Additional File 5, and are repeated without adjustment for these phenomena to allow for assessment of their impact.
We conducted simulations to establish whether such a large effect of fractional location on intensity was plausible. Under our simulations, it is clear that the relative position of the pixel grid over the bead can influence the measured intensity, with the value being greatest when the fractional part of the bead-centre coordinates is close to zero. The range of values that can be achieved is small however, being of the order of 0.2 on the log intensity scale for large bright beads. This value increases as the bead image becomes smaller relative to the pixels (Additional File 6), which would correspond to lower intensity beads, for which such a change might be more important. However this effect is the opposite of the association seen in the real data, suggesting once more that intensity determines location more than location determines intensity.
Ostensibly, this broad result is reassuring. If location were determining intensity to a much larger degree than seen here, then this would be a fundamental problem. Yet even if the majority of the effect is intensity determining location, to allow the loop wherein intensity determines location, which is then allowed to affect intensity, would seem undesirable. The precision of the reported location (to a ten-thousandth of a pixel) would also seem questionable in these circumstances.
For the vast majority of beads the difference between their given and predicted locations is very small, commonly less than one pixel. However, there are often thousands of beads where the deviation is greater than a pixel, and it can commonly be as large as 4 pixels, which is a concern given the beads have a presumed radius of 2-3 pixels and between bead-centre spacings of approximately 6 pixels. Some allowance must be made for irregularities in scanning the array but, given the theoretically regular nature of the grid, this seems to be evidence that beads are being mis-registered, and there is a distinct association between the degree of departure from the grid and residual log-intensity (Additional File 7)
Note that, of the 19 beads, A, D, E, F, G, N, O, P and R map to the wrong position in the red image. Moreover, we see three instances where distinct beads in the green image map to the same bead in the red image: (A,B), (C,D) and (G,H). The consequence of mapping to the wrong bead is clear, with many outliers being created. Bead G should map to the bright bead (in the red channel) to which bead F now maps, but instead maps to a region of low intensity. As a consequence, bead G is taken to exhibit lower intensity than is normal for that bead-type, while bead F exhibits greater intensity than does its cohort. Beads N and O should map to regions of low intensity, but instead return moderately high values as a consequence of their picking out of neighbouring beads. Naturally there are cases, such as P, where the bead-centre maps to the wrong bead in the image, but the two beads exhibit similar intensities so there is little effect.
With localized problems such as this, where frequently only a single bead in a neighbourhood shows a concerning discrepancy between the two channels, it is natural to call upon the redundancy built into the Illumina BeadArray platform and simply discard the beads in question. In these circumstances, correcting a few extra beads would be unlikely to reward adequately the considerable efforts required. The task then is one of identifying these beads for removal, about which we now add a note of caution.
When comparing the coordinates of beads from the two channels, one will immediately be struck by a discrepancy. Since the array is not in exactly the same position in the two images, a shift of one coordinate set will be required, often by a distance of ten or more pixels. However, closer inspection will show that this shift is not constant across the image. The nature of this shift is illustrated in Additional File 8. Two things are clear. First, the array does not occupy the same number of pixels in the red and the green channels, therefore a pixel represents a different physical size between the two channels. Second, for at least one of the channels, the physical size represented by a pixel changes along the array within the channel.
Not only does this have implications for the identification of beads that deviate from the expected grid, but also raises more fundamental questions about the use of a constant number of pixels to calculate foregrounds and backgrounds. Moreover, having seen that the fractional part of a bead's coordinates is influential on the derived intensity, we may ponder what might be the implications of such a changing shift between the two channels.
Sometimes the image registration can go dramatically wrong, resulting in the entire grid of bead-centre locations being positioned incorrectly within the image (Additional File 9). Due to Illumina's random construction method, the consequence of such mis-registration is that the beads have random annotation. The intensity values associated with each bead-type are then a random sample of intensities, and their summarization produces an estimate of the average intensity on the array.
We can automate the re-mapping of bead-types. When the grid of bead IDs is correctly positioned, the variance within each bead-type should be significantly lower than when they are mis-mapped. Being out by a single column or row completely scrambles the annotation, so there should be no improvement until the correct mapping is achieved. We search multiple shifts of the grid and identify when the they are correctly aligned by the drop in variance. Additional File 9 shows, for the mis-registered segment, how the within-bead-type variance drops when the two mis-mapped segments are moved 6 rows vertically. For the two correctly aligned segments, any change in the grid position increases the variance.
We have mentioned already that the decoding process is not perfect and it is common for a small percentage of beads to fail this identification step. The locations of the non-decoded beads fall into two distinct patterns. Many non-decoded beads occur apparently randomly scattered across the array. Others however occur in compact clusters in a manner clearly not independent of one-another, and consistent with a localized technical issue that could have interfered with the decoding step.
Naturally, we must be wary of the possibility that the central bead is not encroaching, but that the satellite beads are showing a signal and that we are simply unable to observed the boundary between them. This may be the case for bead 4, but there are two reasons for believing this not to be the case for the other beads. Firstly, we would expect a local mode of intensity at the satellite bead if this were the case (as there is with bead 4), and secondly we would expect the satellite beads to show similar intensities to beads of a similar type elsewhere on the array.
Figure 6 also shows the difference between the foreground intensity of each of the neighbouring beads and the median intensity for beads of the same type. It is clear that for three of the neighbouring beads (the three over which the fluorescence from the central bead is most visibly encroaching) the foreground intensity is considerably higher than the average for that bead type. Illumina's sharpening step, performed before the foreground intensities are calculated, should go some way to accounting for such encroachment. Given the high residual values for the beads neighbouring a bright bead in Figure 6, which were calculated following this sharpening, it appears that the sharpening may not be adequate for such a purpose (but may still be desirable) . A simple remedy would be to down-weight (or entirely remove) beads that have particularly bright neighbours when summarizing Illumina data. The influence of this bias across an experiment is illustrated in Additional File 10.
Out of the 289 pixels that are considered for the background calculation, the number that are "true background" pixels for a given bead ranges from 117 to 237 for one expression array (4343238066_A_1), with median 149. Note that if we took a 6 × 6 square foreground mask for each bead, more than two-thirds would be considered to have no "true background" pixels within their 289 candidates. However since the outer layer of that 6 × 6 square makes a negative contribution to the foreground calculation, it is reasonable to use the 4 × 4 pixel mask. As anticipated, when there are more true background candidate pixels, the non-robust measure of background tends towards more extreme values. The effect is not large in magnitude however, due to the generally low variability in background values, and so is unlikely to be a concern (see Additional File 11).
When such pixels fall within the background calculation region for a bead, they inevitably end up as one of the five lowest pixels in the region and so contribute to the background score for that bead. The resultant background scores for such beads are lower than they would otherwise be, and the final intensity for the bead is calculated to be higher than it should. There are few such pixels on an array, and their identification and exclusion from calculations would not present a technical challenge, however a simple change to the background calculating rule might present the simplest solution. Currently the mean of the five lowest pixels is used, however using a trimmed mean of a certain number of the lowest pixels or simply using the pixel of a certain rank (equivalent to the median of the n lowest pixels for some value of n) would be alternatives to consider. Figure 7B shows that altering the background calculation to use the median of the five lowest pixels, rather than the mean, decreases the departure from the bead-type type median for the beads in close proximity to the pixel shown in Figure 7A. The influence of this bias across an experiment is illustrated in Additional File 10.
In investigating whether neighbouring beads of the same type can affect the physical chemistry of hybridization to an extent that the registered intensity will also be affected, we find a median of 150 pairs of neighbouring twins (IQR: 142:5 to 159:5) across the 24 sections of the two expression chips. We can detect no significant association between the presence of a neighbour of the same bead-type and observed intensity (Additional File 10), but it should be noted that the range of intensities over which such an effect might be strongest is also that associated with low signal to noise, so we might anticipate having low power to investigate such a phenomenon.
The observed numbers of pairs of such beads far exceed the numbers we expected from theory, or based on a permutation test (Additional File 12). Only two explanations seem plausible for the excessive number of neighbouring twins that we see in the real data. Either the beads are not independently arranged on the array, or the beads are being decoded incorrectly in a manner that leaves them with their neighbour's identity. Should either the physical chemistry (following target depletion) or the various possibilities for explaining the excess neighbour pairs raise concerns, then identification of such bead-pairs is a straightforward task given that the identification of the complete network of neighbouring beads is already a necessary task for various bead-level preprocessing steps. Identification of beads that neighbour similar beads would fall out of this preparation and allow for their down-weighting or exclusion.
In this example the masked bead would normally be included in the summarization. However, its removal results in three additional beads being classed as outliers and excluded from the summarization. This in turn dramatically lowers the summarized intensity and results in the bead type being classed as not expressed, a finding that is corroborated by the replicate sections.
Number of affected beads for expression arrays.
If beads are not providing sound estimates of the quantity that they are designed to measure, yet by chance return an intensity consistent with the true target, then they will have the effect of increasing the perceived number of beads contributing to the estimate and thus will lead to artificially low standard errors and consequently to erroneous results. Since low standard errors are inherently desirable, it is unlikely that concerns will be raised in such cases.
When the intensity is inconsistent with the true target, there is a chance that the standard Illumina analysis will identify these beads as outliers. However it is sub-optimal to use such a broad tool for this purpose, when more targeted options are available. Moreover, we have shown an example where the effect of a miscreant bead is to obscure the outlying status of other beads for which we have no explanation of their abberant intensities. These are the beads for which a general outlier removal step is warranted, and we have shown that a two-pass approach of removing known problems and then applying a general outlier removal can be beneficial.
Summary of results.
Implementation (where implemented)
There is local discordance between the locations in the two channels from two-colour arrays
Between-channel differences in location can be compared to local median differences
If one channel is clearly wrong (from relative grid positions in that channel) then the bead-centre can be remapped from the other channel, else the bead should be dropped
Image cropping out part of the array section, so that values cannot be calculated
Bead-centre coordinates lie without the dimensions of an image, beads apparent on edge of image
Exclude beads with such coordinates
Any text editor can assess the coordinates, while beadarray contains code for the reading and plotting of images on a useful scale. beadarray allows for the masking of beads to be excluded in analyses.
Beads are mostly well-registered, but grid of bead-centres does not align with the image resulting in scrambling of bead-type data
Visual inspection of bead-centres over image. Without access to images, check that the segments have similar extreme x coordinates and that they are equally spaced along the y-axis
In extreme circumstances, bead IDs can be remapped, but usually segments/sections should be excluded
beadarray contains code for the reading and plotting of images, as well as the over-laying of bead-centres. beadarray allows for the masking of beads to be excluded in analyses. Scripts are provided for automatic re-mapping
Neighbouring beads of the same bead-type are a potential concern (albeit unsubstantiated)
Such pairs of beads can be identified and down-weighted or excluded
beadarray contains code for identifying the network neighbourhood. beadarray allows for the masking of beads to be excluded in analyses.
Beads neighbouring clusters of non-decoded beads are more likely to take extreme values
The presence of such clusters can be determined from the .locs file
Exclude/Down-weight beads in a zone about such clusters
beadarray allows for the masking of beads to be excluded in analyses.
Bright beads encroach on neighbours, raising their associated values
Visual inspection of the brightest beads
Bright beads can be identified (by intensity or size using EBImage) and their neighbours down-weighted or excluded
beadarray contains code for the reading and plotting of images, as well as code for identifying neighbours. beadarray allows for the masking of beads to be excluded in analyses.
Abnormally low pixels in the image distort background values and so final intensities
Pixels present with values noticeably lower than the mode
Exclude such pixels, or use a less-sensitive background calculation rule
beadarray contains code to read images, which allow for tests of the distribution of intensities within R. beadarray allows alternative background calculation rules
Multiple bead-centres map to the same location
Text files (or .locs files) can be scanned for neighbours that are unusually close
In two-colour arrays, it may be clear which is the correct bead-centre, else exclude both
Scripts provided to detect departure from predicted bead-centre.
In addition to this, we have demonstrated that in theory the relationship between the position of a bead on the array and the grid of pixels introduced by creating a digital image of the array during scanning can affect the intensity value attributed to that bead. It appears that the majority of the association seen between bead-intensity and within-pixel bead location can be attributed to the intensity of a bead affecting the location at which it is found, yet for a bead to lie away from its anticipated location on the grid we have shown to be indicative of a biased intensity. It may be that different approaches for identifying beads and computing intensities could reduce or eliminate this source of variation.
Finally we have shown that there are a variety of spatial effects that may impact on the reliability of results for individual beads. Particularly bright beads appear to display a 'bleed over' effect in which neighbouring beads are swamped by the signal being emitted from the bright bead. We have also identified isolated pixels with intensities far lower than the modal background for the array. Each of these artificially inflates the intensity attributed to their neighbouring beads, an effect also seen in the neighbourhood of clusters of non-decoded beads.
Yet if they can be identified as suggested in the results and in Table 2, beads can easily be down-weighted or excluded from an analysis. The degree of replication on Illumina arrays often allows for beads to be excluded while still making estimation possible. Unfortunately it is still the case that, on some occasions, entire arrays must be forsaken, so the need for robust experimental design is not negated by these salvage techniques. We have also noted evidence of a lack of independence in the random layout of the BeadArrays. This in itself we have not shown to bias results, however it does suggest that one can not simply trust in an IID random layout of beads to overcome biases introduced by spatial artefacts.
The identification and remedy of the problems we have described are reliant upon access to the raw, bead-level, data. Even if one plans to analyse summarized data using Illumina's software, the routine recording of bead-level data provides a safety net should problems with the arrays be identified or suspected. With such routine access to bead-level data, an automated pipeline for quality assessment could be implemented painlessly, allowing the detection of the artefacts described here, and would provide reassurance regarding an independent summary-level analysis for experimenter, reader, and reviewer alike. Given this we highly recommend the collection and storage of bead-level data.
We thank James Hadfield and Michelle Osbourne for conducting the microarray experiments. We also thank Semyon Kruglyak for useful discussion and advice on Illumina's algorithms, and Matt Ritchie for useful discussions. We acknowledge the support of the University of Cambridge, Cancer Research UK, and Hutchison Whampoa Limited.
This article is published under license to 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.