Method and system for computing and applying a user-defined, global, multi-channel background correction to a feature-based data set obtained from reading a microarray

Abstract

A method and system for estimating a global background-signal correction for each channel of a microarray data set. The method and system of one embodiment of the present invention is directed to a method for calculating background corrected signals for a microarray data set by receiving a non-negative constant and selects a set of low-combined-intensity features from the microarray data set. Based on the low-combined-intensity features, a representation that describes a central-trend of the selected set of low-combined-intensity features is determined in signal-intensity space. The method adjusts the microarray data set parallel to the determined representation based on the non-negative constant.

Description

BACKGROUND OF THE INVENTION

The present invention is related to microarrays. In order to facilitate discussion of the present invention, a general background for particular types of microarrays is provided below. In the following discussion, the terms “microarray,” “microarray,” and “array” are used interchangeably. The terms “microarray” and “microarray” are well known and well understood in the scientific community. As discussed below, a microarray is a precisely manufactured tool which may be used in research, diagnostic testing, or various other analytical techniques to analyze complex solutions of any type of molecule that can be optically or radiometrically detected and that can bind with high specificity to complementary molecules synthesized within, or bound to, discrete features on the surface of a microarray. Because microarrays are widely used for analysis of nucleic acid samples, the following background information on microarrays is introduced in the context of analysis of nucleic acid solutions following a brief background of nucleic acid chemistry.

Deoxyribonucleic acid (“DNA”) and ribonucleic acid (“RNA”) are linear polymers, each synthesized from four different types of subunit molecules. FIG. 1 illustrates a short DNA polymer 100, called an oligomer, composed of the following subunits: (1) deoxy-adenosine 102; (2) deoxy-thymidine 104; (3) deoxy-cytosine 106; and (4) deoxy-guanosine 108. Phosphorylated subunits of DNA and RNA molecules, called “nucleotides,” are linked together through phosphodiester bonds 110-115 to form DNA and RNA polymers. A linear DNA molecule, such as the oligomer shown in FIG. 1, has a 5′ end 118 and a 3′ end 120. A DNA polymer can be chemically characterized by writing, in sequence from the 5′ end to the 3′ end, the single letter abbreviations A, T, C, and G for the nucleotide subunits that together compose the DNA polymer. For example, the oligomer 100 shown in FIG. 1 can be chemically represented as “ATCG.”

The DNA polymers that contain the organization information for living organisms occur in the nuclei of cells in pairs, forming double-stranded DNA helices. One polymer of the pair is laid out in a 5′ to 3′ direction, and the other polymer of the pair is laid out in a 3′ to 5′ direction, or, in other words, the two strands are anti-parallel. The two DNA polymers, or strands, within a double-stranded DNA helix are bound to each other through attractive forces including hydrophobic interactions between stacked purine and pyrimidine bases and hydrogen bonding between purine and pyrimidine bases, the attractive forces emphasized by conformational constraints of DNA polymers. FIGS. 2A-B illustrates the hydrogen bonding between the purine and pyrimidine bases of two anti-parallel DNA strands. AT and GC base pairs, illustrated in FIGS. 2A-B, are known as Watson-Crick (“WC”) base pairs. Two DNA strands linked together by hydrogen bonds forms the familiar helix structure of a double-stranded DNA helix. FIG. 3 illustrates a short section of a DNA double helix 300 comprising a first strand 302 and a second, anti-parallel strand 304.

Double-stranded DNA may be denatured, or converted into single stranded DNA, by changing the ionic strength of the solution containing the double-stranded DNA or by raising the temperature of the solution. Single-stranded DNA polymers may be renatured, or converted back into DNA duplexes, by reversing the denaturing conditions, for example by lowering the temperature of the solution containing complementary single-stranded DNA polymers. During renaturing or hybridization, complementary bases of anti-parallel DNA strands form WC base pairs in a cooperative fashion, leading to reannealing of the DNA duplex.

FIGS. 4-7 illustrate the principle of the microarray-based hybridization assay. A microarray (402 in FIG. 4) comprises a substrate upon which a regular pattern of features is prepared by various manufacturing processes. The microarray 402 in FIG. 4, and in subsequent FIGS. 5-7, has a grid-like 2-dimensional pattern of square features, such as feature 404 shown in the upper left-hand corner of the microarray. Each feature of the microarray contains a large number of identical oligonucleotides covalently bound to the surface of the feature. These bound oligonucleotides are known as probes. In general, chemically distinct probes are bound to the different features of a microarray, so that each feature corresponds to a particular nucleotide sequence.

Once a microarray has been prepared, the microarray may be exposed to a sample solution of target DNA or RNA molecules (410-413 in FIG. 4) labeled with fluorophores, chemiluminescent compounds, or radioactive atoms 415-418. Labeled target DNA or RNA hybridizes through base pairing interactions to the complementary probe DNA, synthesized on the surface of the microarray. FIG. 5 shows a number of such target molecules 502-504 hybridized to complementary probes 505-507, which are in turn bound to the surface of the microarray 402. Targets, such as labeled DNA molecules 508 and 509, that do not contain nucleotide sequences complementary to any of the probes bound to the microarray surface do not hybridize to generate stable duplexes and, as a result, tend to remain in solution. The sample solution is then rinsed from the surface of the microarray, washing away any unbound-labeled DNA molecules. In other embodiments, unlabeled target sample is allowed to hybridize with the microarray first. Typically, such a target sample has been modified with a chemical moiety that will react with a second chemical moiety in subsequent steps. Then, either before or after a wash step, a solution containing the second chemical moiety bound to a label is reacted with the target on the microarray. After washing, the microarray is ready for analysis. Biotin and avidin represent an example of a pair of chemical moieties that can be utilized for such steps.

Finally, as shown in FIG. 6, the bound labeled DNA molecules are detected via optical or radiometric instrumental detection. Optical detection involves exciting labels of bound labeled DNA molecules with electromagnetic radiation of appropriate frequency and detecting fluorescent emissions from the labels, or detecting light emitted from chemiluminescent labels. When radioisotope labels are employed, radiometric detection can be used to detect the signal emitted from the hybridized features. Additional types of signals are also possible, including electrical signals generated by electrical properties of bound target molecules, magnetic properties of bound target molecules, and other such physical properties of bound target molecules that can produce a detectable signal. Optical, radiometric, or other types of instrumental detection produce an analog or digital representation of the microarray as shown in FIG. 7, with features to which labeled target molecules are hybridized similar to 702 optically or digitally differentiated from those features to which no labeled DNA molecules are bound. Features displaying positive signals in the analog or digital representation indicate the presence of DNA molecules with complementary nucleotide sequences in the original sample solution. Moreover, the signal intensity produced by a feature is generally related to the amount of labeled DNA bound to the feature, in turn related to the concentration, in the sample to which the microarray was exposed, of labeled DNA complementary to the oligonucleotide within the feature.

One, two, or more than two data subsets within a data set can be obtained from a single microarray by scanning or reading the microarray for one, two or more than two types of signals. Two or more data subsets can also be obtained by combining data from two different arrays. When optical detection is used to detect fluorescent or chemiluminescent emission from chromophore labels, a first set of signals, or data subset, may be generated by reading the microarray at a first optical wavelength, a second set of signals, or data subset, may be generated by reading the microarray at a second optical wavelength, and additional sets of signals may be generated by detection or reading the microarray at additional optical wavelengths. Different signals may be obtained from a microarray by radiometric detection of radioactive emissions at one, two, or more than two different energy levels. Target molecules may be labeled with either a first chromophore that emits light at a first wavelength, or a second chromophore that emits light at a second wavelength. Following hybridization, the microarray can be read at the first wavelength to detect target molecules, labeled with the first chromophore, hybridized to features of the microarray, and can then be read at the second wavelength to detect target molecules, labeled with the second chromophore, hybridized to the features of the microarray. In one common microarray system, the first chromophore emits light at a near infrared wavelength, and the second chromophore emits light at a yellow visible-light wavelength, although these two chromophores, and corresponding signals, are referred to as “red” and “green.” The data set obtained from reading the microarray at the red wavelength is referred to as the “red signal,” and the data set obtained from reading the microarray at the green wavelength is referred to as the “green signal.”While it is common to use one or two different chromophores, it is possible to use one, three, four, or more than four different chromophores and to read a microarray at one, three, four, or more than four wavelengths to produce one, three, four, or more than four data sets. With the use of quantum-dot dye particles, the emission is tunable by suitable engineering of the quantum-dot dye particles, and a fairly large set of such quantum-dot dye particles can be excited with a single-color, single-laser-based excitation.

FIG. 8 shows a small region of a scanned image of a microarray containing an image of a single feature. In FIG. 8, the small region of the scanned image comprises a grid, or matrix, of pixels, such as pixel 802. In FIG. 8, the magnitude of the signal scanned from the small region of the surface of a microarray spatially corresponding to a particular pixel in the scanned image is indicated by a kind of gray scaling. Pixels corresponding to high-intensity signals, such as pixel 804, are darkly colored, while pixels having very low signal intensities, such as pixel 802, are not colored. The range of intermediate signal intensities is represented, in FIG. 8, by a generally decreasing density of crosshatch lines within a pixel. In FIG. 8, there is a generally disc-shaped region in the center of the region of the scanned image of the microarray that contains a high proportion of high-intensity pixels. Outside of this central, disc-shaped region corresponding to a feature, the intensities of the pixels fall off relatively quickly, although pixels with intermediate intensities are found, infrequently, even toward the edges of the region of the scanned image, relatively distant from the obvious central, disc-shaped region of high-intensity pixels that corresponds to the feature.

In general, data sets collected from microarrays comprise an indexed set of numerical signal intensities associated with pixels. The pixel intensities range over the possible values for the size of the memory-storage unit employed to store the pixel intensities. In many current systems, a 16-bit word is employed to store the intensity value associated with each pixel, and a data set can be considered to be a 2-dimensional array of pixel-intensity values corresponding to the 2-dimensional array of pixels that together compose a scanned image of a microarray.

FIG. 9 shows a 2-dimensional array of pixel-intensity values corresponding to a portion of the central, disc-shaped region corresponding to a feature in the region of a scanned image of a microarray shown on FIG. 8. In FIG. 9, for example, pixel intensity 902 corresponds to pixel 806 in FIG. 8.

Features on the surface of a microarray may have various different shapes and sizes, depending on the manufacturing process by which the microarray is produced. In one important class of microarrays, features are tiny, disc-shaped regions on the surface of the microarray produced by ink-jet-based application of probe molecules, or probe-molecular-precursors, to the surface of the microarray substrate. FIG. 10 shows an idealized representation of a feature, such as the feature shown in FIG. 8, on a small section of the surface of a microarray. FIG. 11 shows a graph of pixel-intensity values versus position along a line bisecting a feature in the scanned image of the feature. For example, the graph shown in FIG. 11 may be obtained by plotting the intensity values associated with pixels along lines 1002 or 1004 in FIG. 10. Consider a traversal of the pixels along line 1002 starting from point 1006 and ending at point 1008. In FIG. 11, points 1106 and 1108 along the horizontal axis correspond to positions 1006 and 1008 along line 1002 in FIG. 10. Initially, at positions well removed from the central, disc-shaped region of the feature in 1010, the scanned signal intensity is relatively low. As the central, disc-shaped region of the feature is approached, along line 1002, the pixels intensities remain at a fairly constant, background level up to point 1012, corresponding to point 1112 in FIG. 11. Between points 1012 and 1014, corresponding to points 1112 and 1114 in FIG. 11, the average intensity of pixels rapidly increases to a relatively high intensity level 1115 at a point 1014 coincident with the outer edge of the central, disc-shaped region of the feature. The intensity remains relatively high over the central, disc-shaped region of the feature 1116, and begins to fall off starting at point 1018, corresponding to point 1118 in FIG. 11, at the far side of the central, disc-shaped region of the feature. The intensity rapidly falls off with increasing distance from the central, disc-shaped region of the feature until again reaching a relatively constant, background level at point 1008, corresponding to point 1108 in FIG. 11. The exact shape of the pixel-intensity-versus-position graph, and the size and shape of the feature, are dependent on the particular type of microarray and molecular-array substrate, chromophore or a radioisotope used to label target molecules, experimental conditions to which the microarray is subjected, the molecular-array reader used to scan a microarray, and on data processing components of the molecular-array reader and an associated computer that produce the scanned image and pixel-intensity data sets. For example, with some type of array manufacture processes or with different hybridization and washing protocols, the features may resemble donuts, or even more irregular blobs.

The background signal generated during scanning regions of the surface of a microarray outside of the areas corresponding to features arises from many different sources, including contamination of the molecular-array surface by fluorescent or radioactively labeled or naturally radioactive compounds, fluorescence or radiation emission from the molecular-array substrate, dark signal generated by the photo detectors in the molecular-array reader, and many other sources. When this background signal is measured on the portion of the array that is outside of the areas corresponding to a feature, it is often referred to as the “local” background signal.

An important part of molecular-array data processing is a determination of the background signal that needs to be subtracted from a feature. With appropriate background-subtraction, it is possible to distinguish low-signal features from no-signal features and to calculate accurate and reproducible log ratios between multi-channel and/or inter-array data. The sources of background signal that appear in the local background region may be identical to the sources of background signal that occur on the feature itself; that is, the signal represented in the local background region may be correlated with the signal that arises from the specific labeled target hybridized to probes on that feature. In this case, it is appropriate to use the signal from the local background region as the best estimate of the background to subtract from that feature.

FIG. 12 illustrates a currently employed technique for measuring the local background signal for a feature. FIG. 12 corresponds to the small region of the scanned image of a microarray shown on FIG. 8. Initial pixel-based coordinates for the center of the feature can be estimated from manufacturing data for the microarray and from a number of scanned-image processing techniques. Using these initial pixel-based coordinates for the center of the feature, the integrated intensities of disc-shaped regions with increasing radii centered at those coordinates can be computed to determine, by a decrease in integrated intensities, the outer edge 1202 of the central, disc-shaped feature region. An intermediate region, in which the integrated pixel intensities rapidly fall off with increasing radius, corresponding to the regions in FIG. 11 between points 112 and 114 and between points 1118 and 1108, can be determined to provide the outer boundary 1204 of a region of interest (“ROI”) surrounding and including the central, disc-shaped region of the feature. Finally, an annulus lying between the outer edge of the ROI 1204 and a somewhat arbitrary outer background circumference 1206 is considered to be the background region for the feature, and the integrated intensity of this background region 1208, divided by the area of the background region, is taken to be the background signal for the entire region comprising the feature region, feature ROI, and background annulus. Alternatively, the locations and sizes of feature regions may be known in advance of the image processing stage, based on array manufacturing data and other information, and so the ROI may not need to be determined by a method such as the method described above. Thus, a current technique for background signal estimation is based on a local method involving determining an integrated signal intensity for an annulus surrounding the ROI disc associated with a feature, and determining a background-signal intensity per image area. The estimated local background signal for a feature is the background-signal intensity per image area, and is subtracted from the normalized raw feature signal, to produce a background-subtracted feature signal. A feature-based data set includes background-subtracted data subset, for each signal scanned, comprising feature signals or raw feature signals.

In general, for a large class of molecular-array-based experiments, the logarithm of the ratio of two feature signals, measured in two channels, for each feature, referred to as the log ratio, is calculated as a measure of the differential concentrations of target molecules in those two sample solutions that bind to the feature. A plurality of channels may include both intra-array channels, two or more signal channels scanned from a single array, and may also include inter-array channels, one or more signal channels scanned from two or more arrays. Unfortunately, when both signals are relatively weak, or of low magnitude, the log ratio can be extremely sensitive to slight perturbations in the relative weak intensities of the two signals. Imprecise background signal correction can lead to anomalous log ratios and spurious results based on the anomalous log ratios. Manufacturers and designers of microarrays and microarray readers, as well as researchers and diagnosticians who use microarrays in experimental and commercial settings, have recognized the need for accurate background subtraction methods, in order to obtain more accurate and reproducible gene expression data.

SUMMARY OF THE INVENTION

Various embodiments of the present invention are directed toward estimating and correcting background signals present in microarray data. One embodiment of the present invention is directed to a method and system for calculating background corrected signals for a microarray data set by receiving a non-negative constant and selects a set of low-combined-intensity features from the microarray data set. Based on the low-combined-intensity features, a representation that describes a central-trend of the selected set of low-combined-intensity features is determined. The method adjusts the microarray data set parallel to the determined representation based on the non-negative constant.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a short DNA polymer.

FIGS. 2A-B illustrate the hydrogen bonding between the purine and pyrimidine bases of two anti-parallel DNA strands.

FIG. 3 illustrates a short section of a DNA double helix comprising a first strand and a second, anti-parallel strand.

FIG. 4 illustrates a grid-like, two-dimensional pattern of square features.

FIG. 5 shows a number of target molecules hybridized to complementary probes, which are in turn bound to the surface of the microarray.

FIG. 6 illustrates the bound labeled DNA molecules detected via optical or radiometric scanning.

FIG. 7 illustrates optical, radiometric, or other types of scanning produced by an analog or digital representation of the microarray.

FIG. 8 shows a small region of a scanned image of a microarray containing the image of a single feature.

FIG. 9 shows a 2-dimensional array of pixel-intensity values corresponding to a portion of the central, disc-shaped region corresponding to a feature in the region of a scanned image of a microarray shown on FIG. 8.

FIG. 10 shows an idealized representation of a feature, such as the feature shown in FIG. 8, on a small section of the surface of a microarray.

FIG. 11 shows a graph of pixel-intensity values versus position along a line bisecting a feature in the scanned image of the feature.

FIG. 12 illustrates a currently employed technique for measuring the background signal for a feature.

FIG. 13 illustrates a currently available approach to residual, global background signal correction of a two-channel, feature-based data set.

FIGS. 14-15 illustrate low-intensity anomalies that occur in graphically analyzed data sets resulting from incorrect or imprecise global background-signal correction by the technique shown in FIG. 13.

FIG. 16 illustrates a better approach to correcting residual background signals.

FIG. 17 illustrates an over-corrected data set resulting from local background-signal subtraction.

FIG. 18 is a high-level flow diagram describing the method of one embodiment of the present invention.

FIGS. 19A-B show a feature having uniform intensity distribution and a feature have non-uniform intensity distribution.

FIG. 20 illustrates the central trend within a data set.

FIGS. 21A-B, 22 A-B, and 23-24 illustrate a rank-order method for selecting features from a set of filtered features.

FIG. 25 illustrates the fitting of a line to the low-combined-intensity central-trend features.

FIG. 26 illustrates the principle of the MAD technique.

FIGS. 27-28 illustrate two methods that can be employed to filter the low-combined-intensity, central-trend features with respect to the fitted line.

FIGS. 29-30 illustrate the construction of an ideal background data point.

FIG. 31 illustrates a final refinement technique for calculating the ideal characteristic background feature.

FIG. 32 is a control-flow diagram of the routine “Correct C1 and C2 channel data.”

FIGS. 33A-B illustrate moving a C1/C2 data set having a best-fit line slope greater than “1.”

FIGS. 34A-B illustrate moving a C1/C2 data set having a best-fit line slope less than “1.”

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

One embodiment on the present invention is directed to a method and system for determining global-background-signal-based background corrections for each channel in a feature-based, molecular-array data set. The term global refers to an identical correction applied to all features regardless of position on the microarray or signal level. In one embodiment, described below, the method provides separate global, residual-background-signal corrections for each of two channels C1 and C2. This method is useful for current molecular-array-based experiments in which two different chromophores, radioactive labels, or other types of labels are employed and scanned in two different channels in a molecular reader. However, the present invention is not restricted to determining global, residual-background-signal corrections for two channels, but can be extended to additional channels when more than two labels are employed in molecular-array-based experiments and detected by a microarray reader. The present invention is not restricted to correcting data from a single array, or intra-array data. There are two broad categories of potential corrections for data obtained from multiple arrays, or inter-array data. The first category involves two arrays that are each single-channel arrays, and a need to correctly compare C1 signal from a first array with C1 signal from a second array. The second category involves multi-channel arrays, and a need, for example, to compare C1 signal from a first array with C2 signal from a second array. The present invention is not restricted to correcting data from background-subtracted signals. An embodiment of the method of the present invention can be applied to raw signals, before any background-subtraction is done, and thus includes both background-signal subtraction and a global, residual background correction method in a single step.

The described embodiment of the present invention is directed to identifying C1 and C2 residual background-signal corrections. When the C1-signal and C 2-signal intensities of each feature are plotted in a C1/C2 plot, the data points corresponding to features are generally distributed about a central-trend line or curve. In general, there is an apparent cluster of smallest-intensity features at a low-intensity region of the distribution, or lowest-intensity point of the central-trend line or curve. In the described embodiment, the C1 and C2 residual bacground-signal corrections are equivalent to the x′ and y′ coordinates for an ideal, characteristic, low intensity data point within the cluster of smallest-intensity features in a C1/C2 plot. The C1 and C2 feature intensities for the features in the data set can be corrected to translate the ideal, characteristic, low intensity data point to the origin of the C1/C2 plot, or equivalently, the magnitudes of the C1 and C2 residual background-signal corrections correspond to the x′ and y′ coordinates for the ideal, characteristic, low intensity data point in a C1/C2 plot. One embodiment of the present invention is described in four subsections that follow: (1) additional information about microarrays; (2) an overview of the method of one embodiment of the present invention, presented with reference to FIGS. 13-17; (3) a flow-control-diagram and illustration-based description of the method of the embodiment of the present invention; and (4) a Matlab-like pseudocode implementation of the embodiment.

Additional Information About Microarrays

A microarray may include any one-, two-or three-dimensional arrangement of addressable regions, or features, each bearing a particular chemical moiety or moieties, such as biopolymers, associated with that region. Any given microarray substrate may carry one, two, or four or more microarrays disposed on a front surface of the substrate. Depending upon the use, any or all of the microarrays may be the same or different from one another and each may contain multiple spots or features. A typical microarray may contain more than ten, more than one hundred, more than one thousand, more ten thousand features, or even more than one hundred thousand features, in an area of less than 20 cm² or even less than 10 cm². For example, square features may have widths, or round feature may have diameters, in the range from a 10 μm to 1.0 cm. In other embodiments each feature may have a width or diameter in the range of 1.0 μm to 1.0 mm, usually 5.0 μm to 500 μm, and more usually 10 μM to 200 μm. Features other than round or square may have area ranges equivalent to that of circular features with the foregoing diameter ranges. At least some, or all, of the features may be of different compositions (for example, when any repeats of each feature composition are excluded the remaining features may account for at least 5%, 10%, or 20% of the total number of features). Inter-feature areas are typically, but not necessarily, present. Inter-feature areas generally do not carry probe molecules. Such inter-feature areas typically are present where the microarrays are formed by processes involving drop deposition of reagents, but may not be present when, for example, photolithographic microarray fabrication processes are used. When present, interfeature areas can be of various sizes and configurations.

Each microarray may cover an area of less than 100 cm², or even less than 50 cm², 10 cm² or 1 cm². In many embodiments, the substrate carrying the one or more microarrays will be shaped generally as a rectangular solid having a length of more than 4 mm and less than 1 m, usually more than 4 mm and less than 600 mm, more usually less than 400 mm; a width of more than 4 mm and less than 1 m, usually less than 500 mm and more usually less than 400 mm; and a thickness of more than 0.01 mm and less than 5.0 mm, usually more than 0.1 mm and less than 2 mm and more usually more than 0.2 and less than 1 mm. Other shapes are possible, as well. With microarrays that are read by detecting fluorescence, the substrate may be of a material that emits low fluorescence upon illumination with the excitation light. Additionally in this situation, the substrate may be relatively transparent to reduce the absorption of the incident illuminating laser light and subsequent heating if the focused laser beam travels too slowly over a region. For example, a substrate may transmit at least 20%, or 50% (or even at least 70%, 90%, or 95%), of the illuminating light incident on the front as may be measured across the entire integrated spectrum of such illuminating light or alternatively at 532 nm or 633 nm.

Microarrays can be fabricated using drop deposition from pulsejets of either polynucleotide precursor units (such as monomers) in the case of in situ fabrication, or the previously obtained polynucleotide. Such methods are described in detail in, for example, U.S. Pat. No. 6,242,266, U.S. Pat. No. 6,232,072, U.S. Pat. No. 6,180,351, U.S. Pat. No. 6,171,797, U.S. Pat. No. 6,323,043, U.S. patent application Ser. No. 09/302,898 filed Apr. 30, 1999 by Caren et al., and the references cited therein. Other drop deposition methods can be used for fabrication, as previously described herein. Also, instead of drop deposition methods, photolithographic microarray fabrication methods may be used. Interfeature areas need not be present particularly when the microarrays are made by photolithographic methods as described in those patents.

A microarray is typically exposed to a sample including labeled target molecules, or, as mentioned above, to a sample including unlabeled target molecules followed by exposure to labeled molecules that bind to unlabeled target molecules bound to the microarray, and the microarray is then read. Reading of the microarray may be accomplished by illuminating the microarray and reading the location and intensity of resulting fluorescence at multiple regions on each feature of the microarray. For example, a scanner may be used for this purpose, which is similar to the AGILENT MICROARRAY SCANNER manufactured by Agilent Technologies, Palo Alto, Calif. Other suitable apparatus and methods are described in published U.S. patent applications 20030160183A1, 20020160369A1, 20040023224A1, and 20040021055A, as well as U.S. Pat. No. 6,406,849. However, microarrays may be read by any other method or apparatus than the foregoing, with other reading methods including other optical techniques, such as detecting chemiluminescent or electroluminescent labels, or electrical techniques, for where each feature is provided with an electrode to detect hybridization at that feature in a manner disclosed in U.S. Pat. No. 6,251,685, and elsewhere.

A result obtained from reading a microarray, followed by application of a method of the present invention, may be used in that form or may be further processed to generate a result such as that obtained by forming conclusions based on the pattern read from the microarray, such as whether or not a particular target sequence may have been present in the sample, or whether or not a pattern indicates a particular condition of an organism from which the sample came. A result of the reading, whether further processed or not, may be forwarded, such as by communication, to a remote location if desired, and received there for further use, such as for further processing. When one item is indicated as being remote from another, this is referenced that the two items are at least in different buildings, and may be at least one mile, ten miles, or at least one hundred miles apart. Communicating information references transmitting the data representing that information as electrical signals over a suitable communication channel, for example, over a private or public network. Forwarding an item refers to any means of getting the item from one location to the next, whether by physically tran-sporting that item or, in the case of data, physically transporting a medium carrying the data or communicating the data.

As pointed out above, microarray-based assays can involve other types of biopolymers, synthetic polymers, and other types of chemical entities. A biopolymer is a polymer of one or more types of repeating units. Biopolymers are typically found in biological systems and particularly include polysaccharides, peptides, and polynucleotides, as well as their analogs such as those compounds composed of, or containing, amino acid analogs or non-amino-acid groups, or nucleotide analogs or non-nucleotide groups. This includes polynucleotides in which the conventional backbone has been replaced with a non-naturally occurring or synthetic backbone, and nucleic acids, or synthetic or naturally occurring nucleic-acid analogs, in which one or more of the conventional bases has been replaced with a natural or synthetic group capable of participating in Watson-Crick-type hydrogen bonding interactions. Polynucleotides include single or multiple-stranded configurations, where one or more of the strands may or may not be completely aligned with another. For example, a biopolymer includes DNA, RNA, oligonucleotides, and PNA and other polynucleotides as described in U.S. Pat. No. 5,948,902 and references cited therein, regardless of the source. An oligonucleotide is a nucleotide multimer of about 10 to 100 nucleotides in length, while a polynucleotide includes a nucleotide multimer having any number of nucleotides.

As an example of a non-nucleic-acid-based microarray, protein antibodies may be attached to features of the microarray that would bind to soluble labeled antigens in a sample solution. Many other types of chemical assays may be facilitated by microarray technologies. For example, polysaccharides, glycoproteins, synthetic copolymers, including block copolymers, biopolymer-like polymers with synthetic or derivitized monomers or monomer linkages, and many other types of chemical or biochemical entities may serve as probe and target molecules for microarray-based analysis. A fundamental principle upon which microarrays are based is that of specific recognition, by probe molecules affixed to the microarray, of target molecules, whether by sequence-mediated binding affinities, binding affinities based on conformational or topological properties of probe and target molecules, or binding affinities based on spatial distribution of electrical charge on the surfaces of target and probe molecules.

Scanning of a microarray by an optical scanning device or radiometric scanning device generally produces an image comprising a rectilinear grid of pixels, with each pixel having a corresponding signal intensity. These signal intensities are processed by a microarray-data-processing program that analyzes data scanned from an microarray to produce experimental or diagnostic results which are stored in a computer-readable medium, transferred to an intercommunicating entity via electronic signals, printed in a human-readable format, or otherwise made available for further use. Microarray experiments can indicate precise gene-expression responses of organisms to drugs, other chemical and biological substances, environmental factors, and other effects. Microarray experiments can also be used to diagnose disease, for gene sequencing, and for analytical chemistry. Processing of microarray data can produce detailed chemical and biological analyses, disease diagnoses, and other information that can be stored in a computer-readable medium, transferred to an intercommunicating entity via electronic signals, printed in a human-readable format, or otherwise made available for further use.

Overview

In some cases, it is more appropriate to use a global method to determine the background signal and then to subtract that global background signal from all features of an array. One example is to calculate an average or median of all the local background regions on an array. Another method is to determine a minimum feature or local-background-region signal. Another method is to use a set of features, referred to as “negative controls,” that are designed to not bind labeled target in a specific manner.

To determine the correctness of a background subtraction method, it is common to analyze data from self-self arrays. A self-self experiment generally involves taking one sample and splitting it for labeling. A first portion of the sample is labeled with a first label and an equal, second portion of the sample is labeled with a second label. The two portions are then applied to a single array for hybridization, washing, and scanning. Three types of graphical analysis of such an array should yield the following three results, respectively, after correct background subtraction: (1) a C1/C2 plot plotted linearly in both dimensions is generally approximately linear and symmetric; (2) a C1/C2 plot plotted in log scale is generally approximately linear and symmetric; and (3) a plot of the log of the ratio C1/C2 versus C1 or C2 should be approximately horizontal. The slope of the lines from first and second graphical analyses, described above, and the y-intercept from the third graphical analysis, described above, reflect a constant that is needed for dye-normalization. Such a constant reflects different efficiencies in labeling for the C1 and C2 labels and different efficiencies in scanner optics. Once the data is dye-normalized, the data generally yields a slope of 1.0 for the C1/C2 plots and ay intercept of 0 for the plot of the log of the ratio C1/C2 versus C1 or C2. For the purpose of this discussion, the dye-normalization constant can be applied to features in a signal-independent fashion. In many observed cases, dye-normalization is often signal-dependent and curve-fitting may be needed to achieve optimal results. Indeed, for cases where a curve-fitting algorithm is employed for dye-normalization, the results will be more robust if correct background subtraction has first been applied. When background-signal correction is incorrect or imprecise, the ideal, expected results from the above-described graphical analyses, and expected results for similar analytic approaches for the curve-fitting techniques needed for non-linear cases, are not observed. Instead, the graphical analyses produce anomalous trends and features.

FIG. 13 illustrates one approach to residual, global background signal correction of a two-channel, feature-based data set. FIG. 13 is a C1/C2 plot of the features in a hypothetical data set. In FIG. 13, the horizontal x axis 1304 corresponds to the intensity of a feature in the C2 channel, and the vertical y axis 1306 corresponds to the feature intensity in the C1 channel. In FIG. 13, a data point, such as data point 1302, is plotted with xy coordinates equal to a corresponding feature's C 2-signal and Cl-signal intensities. The data points corresponding to features often form a cloud-like distribution about a line or curve fit, by well-known mathematical curve-fitting techniques, to the distribution. For example, in FIG. 13, a line 1308 has been fit to the distribution of points. C2 and C1 may represent feature intensities measured at green and red wavelengths, respectively; or may represent feature intensities measured at the same wavelength from two different arrays, or feature intensities measured for two different labels via radiometric, magnetic, or electrical scanning techniques

In general, the feature intensities within a data set span a relatively large range of values, from very high intensities down to a lowest set of intensities generally somewhat above the 0-level intensity, for non-background-subtracted signals. In the case of background-subtracted signals, the lowest set of signals may have negative values. In a self-self experiment, the data points should fall symmetrically about a line with a slope that corresponds to a constant C1/C2 ratio that reflects the relative target labeling efficiencies and relative efficiencies by which the C1 and C2 signals are measured by a microarray reader. Again, as noted above, the number of C1 and C2 labeled target molecules resident within any particular feature in a microarray produced in a self-self experiment should be nearly identical. Even in a differential expression experiment, the distribution of data points generally exhibits a central trend, to which a line or curve can be mathematically fit.

In one currently employed background-signal correction technique, the best-fit line, 1308 in FIG. 13, is extended until it intersects with the x axis or y axis. For example, in FIG. 13, the best-fit line is extended downward and leftward to the point 1310, the x-axis intercept. The data subset corresponding to the axis intercepted, in the case of FIG. 13, the x axis corresponding to the C2 data subset, is then corrected by applying the magnitude of the x-axis correction 1312 to the data subset to translate the intersection point 1310 to the origin of the C1/C2 plot. In general, the sign of the x or y coordinate of the intercept is reversed, and the x or y coordinate added to the feature signals. However, this technique is inadequate. In general, residual background signals are present in both C1 and C2 data subsets, so that the relative position of the distribution of data points is offset in both x and y directions with respect to the coordinate system. Moreover, extending the central-trend curve to an x-axis or y-axis intersection is an arbitrary construct, without sound theoretical basis.

FIGS. 14-15 illustrate low-intensity anomalies that occur in graphically analyzed data sets resulting from incorrect or imprecise global background-signal correction by the technique shown in FIG. 13. FIG. 14 shows a representation of a log-ratio-versus-C 2-signal intensity plot for features of a data set. Note that, in FIG. 14, and in many subsequent figures, a curve is shown to represent the central trend of a distribution of data points that would be plotted from an actual data set. The data set plotted in FIG. 14 would produce a C1/C2 plot similar to that shown in FIG. 13, with a C1/C2 best-fit line having a slope less than 1. A log-ratio-versus-C2-signal intensity plot would be expected to produce an essentially horizontal, or constant, line, offset by the value log(m) (1402 in FIG. 14), where m is the slope of the C1/C2 best-fit line in a C1/C2 plot, from the horizontal axis 1404. For large signal intensity ratios, the log-ratio-versus-C 2-signal intensity plot does approximate a horizontal line with the expected offset from the x axis, as at data point 1406. However, as the combined-signal feature intensity decreases, toward the y axis, the log-ratio-versus-C2-signal intensity central-trend curve begins to dramatically curve upward, crossing the x axis at point 1410, and asymptotically approaching the y axis for very small values. Plots of log(C1) versus log (C2 ) exhibit similar, low-intensity anomalies. As shown in FIG. 15, a log(C1)-versus-log-(C2 ) is also expected to show a linear central-trend line for a data set with a relatively constant C1/C2 ratio. But, for low intensity features, the central-trend line 1502, rather than continuing 1504 at a constant slope towards an x-axis or y-axis intercept, is often observed to curve away from line 1502 in either a positive direction 1506 or negative direction 1508.

An important contribution to the anomalous, low-intensity departure of observed central-trend curves from expected central-trend curves for log-ratio-versus-C2-signal and log(C1)-versus-log-(C2 ) plots, discussed above with reference to FIGS. 14 and 15, is the flawed approach to global, residual background correction, noted above with reference to FIG. 13. If, after global, residual background correction, residual background signals remain within the data set, and if the residual background signals are of different magnitudes for the two channels, then, for low intensity features, as the feature intensity in one channel approaches the corrected 0-signal level, the feature intensity in the other channel may approach the value of the difference between the residual background signals in the two channels. The log ratio then begins to rapidly depart from a general, constant value, veering non-linearly either towards plus infinity or towards minus infinity. In general, the upwards curving seen in FIG. 14 and FIG. 15 occurs when C1 is under-corrected for background signal, relative to C2, and downwards curving, as seen in FIG. 15, occurs when C2 is under-corrected for background signal, relative to C1.

FIG. 16 illustrates a better approach to correcting residual background signals. In FIG. 16, a central-trend line 1602 is fit to a distribution of points, such as point 1604, in a C1/C2 plot. However, in the better approach, instead of arbitrarily extending the central-trend line towards an x-axis or y-axis intercept, the location of an ideal, characteristic, low intensity data point 1606 is estimated based on the lowest intensity, non-outlying data points clustered about the central-trend line. The x and y coordinates 1608 and 1610 of this ideal, characteristic, low intensity data point represents C2 and C1 global residual background-signal corrections. Application of the C2 and C1 global residual background-signal corrections has the effect of moving the position of the ideal, characteristic, low intensity data point to the origin of the C1/C2 plot. Unlike the technique described with respect to FIG. 13, in which one data subset is shifted relative to the other data subset, the better approach illustrated in FIG. 16 involves calculating both C2 and C1 global residual background-signal corrections, based upon perpendicular projections to the axes rather than a linear extrapolation, and applying both C2 and C1 global residual background-signal corrections to their respective data subsets. The location of the ideal, characteristic, low intensity data point is a good approximation for a data point corresponding to an ideal, no-background-signal-in-either-channel feature, and its relative position with respect to the x and y axes reflects the respective global, residual background signals in the C2 and C1 data subsets. By contrast, the x-axis intersection 1609 of the central-trend line is reflective of neither the C2 nor C1 global residual background-signals.

It should be noted that global, residual background signals in the C2 and C1 data subsets may be positive or negative for local-background-subtracted data sets. Local-background subtraction may lead to over-correction in one or both channels, in which case positive global residual background-signal corrections may need to be applied in one or both channels. For a “positive” global residual background-signal correction, adding back the residual value is needed, while for a “negative” global background-signal correction, subtracting the residual value is needed. FIG. 17 illustrates an over-corrected data set resulting from local background-signal subtraction. In FIG. 17, the ideal, characteristic, low intensity data point 1702 falls within the fourth quadrant, below the x axis. In this case, a positive C1 global, residual background-signal correction and a negative C2 global residual background-signal correction need to be applied to the data set in order to bring the ideal, characteristic, low intensity data point 1702 to the origin of the C1/C2 plot. Had the location of the ideal, characteristic, low intensity data point fallen in the third quadrant 1704, then positive C1 and C2 global residual background-signal corrections would need to be applied to the data set.

Note that, in quasi-self-self experiments, in which multiple 1-channel arrays are exposed to the same sample solution and scanned, C1/Cl′ plots can be employed to identify the position of an ideal, characteristic, low intensity data point to allow for calculation of separate, global, residual background-signal corrections for the C1 and C1′ data subsets. In both self-self, and quasi-self-self experiments, the distribution of data points about a central-trend line tends to be symmetrical. In differential experiments, although not fully symmetrical, the distribution nonetheless clusters about a central-trend line or curve and the subject of this invention will provide an effective background correction method. A quasi-differential expression experiment can result from comparing a signal set from one array with one sample type versus a second array hybridized with a second sample type. This can occur with either single or multi-channel array data. The subject of this invention will also provide an effective background correction method for these cases.

Flow-Control Diagram and Graphical Description of One Embodiment

FIG. 18 is a high-level flow diagram describing the method of one embodiment of the present invention. This method will be described with reference to FIGS. 18 and concurrent reference to subsequent figures that provide greater detail for many of the steps shown in FIG. 18.

In the first step 1802 of the two-channel background correction method, a two-channel, feature-based molecular-array data set is received. As discussed earlier, this data set can comprise two data subsets corresponding to different channels of one array, or this data set can comprise two data subsets corresponding to signals obtained from different arrays. The feature intensities of the data set can be either raw intensities or local-background-subtracted intensities. Next, in the step 1804, those features in the data set not marked as control features are selected. Initially, control features are filtered out, or rejected, because they may not to exhibit C1/C2 ratios close to the central trend for C1/C2 ratios within the data set.

In step 1806, the selected non-control features are filtered to remove features with a saturation level above a threshold saturation level and features with non-uniform pixel-intensity distributions. As discussed above with reference to FIG. 9, the scanned image of a microarray consists of a 2-dimensional array of pixel-intensity values, commonly stored in 16-bit words, and therefore ranging from 0 to 65535. A pixel having an intensity value of 65535 is considered to be saturated, because all measured intensity values greater that 65335 are encoded as the maximum value 65535. When more than a threshold percentage of the pixels within the area corresponding to a feature are saturated, then the feature is considered to be saturated. In other words, the true intensity of the feature is not reflected in the intensity value integrated over the pixels within the feature area. In one embodiment of the present invention, the saturation-level threshold is 5% of the feature pixels, and features having more than 5% saturated pixels are removed from the set.

A second filter, carried out in step 1806, is designed to remove features with non-uniform intensity distributions over the area of the feature. FIGS. 19A-B show a feature having uniform intensity distribution and a feature have non-uniform intensity distribution. In FIG. 19A, the disc-shaped area of the image of a feature 1902 is shown with cross-hatching indicating a uniform distribution of moderate intensity values over the entire area of the image of the feature. By contrast, in FIG. 19B, a crescent-shaped portion 1904 of the area of a feature has medium intensity values, while the remaining portion of a feature, 1906, has low intensity values. The signal intensities within the feature shown in FIG. 19B are non-uniformly distributed. Non-uniform distribution of intensities can be detected in a number of different of ways. The statistical variance of pixel intensities within the area of image of a feature can be computed, and features with pixel-intensity variances greater than a threshold variance can be considered to have non-uniform pixel-intensity distributions. The threshold may be a function of biological or electronic noise, or noise on the microarray due to chemical or hybridization processes, for example. Alternatively, as shown in FIG. 19B, a centroid for the pixel intensity distribution 1908 can be computed, and the location of the centroid compared to the location of the geometric center 1910 of the image of the feature. When the distance 1912 between the centroid of pixel-intensity and the geometric center is greater than a threshold distance value, the feature can be considered to have a non-uniform pixel-intensity distribution. Alternatively, both the mean and the median signal statistics can be computed and a feature can be considered to have a non-uniform pixel-intensity distribution if the difference between the mean signal and the median signal exceeds a threshold or if the difference divided by either the mean or the median signal exceeds a threshold. Many other alternative techniques can be employed in order to classify features having uniform and non-uniform pixel-intensity distributions. In step 1806, features having non-uniform pixel distributions are removed from the selected set of non-control features. As noted above, features in certain types of arrays may not be disc-shaped, and techniques used to compute a metric of uniformity may need to be tailored specifically to the particular feature shapes present in the arrays.

In step 1808, features that fall within a central trend within the data set are chosen from the set of filtered features produced in step 1806. Such features exhibit generally an equal, negligibly different, relation between the signals in the two channels.

FIG. 20 illustrates the central trend within a data set. In FIG. 20, as in FIGS. 13, 16, and 17, described above, each feature is plotted with respect to the features′ signal-intensities in the C2 and C1 channels. Features close to the line 2002 are features having a C1/C2 ratio close to a characteristic C1/C2 ratio for the data set. Features further from the line, such as feature 2004, have C1/C2 ratios that markedly depart from the characteristic C1/C2 ratio for the data set. In differential experiments, a certain portion of the data points are expected to fall outside of close proximity to the central-trend, reflecting different signals arising from different concentrations of target molecules in the two sample solutions, and the overall distribution of data points may be somewhat asymmetrical. However, even in arrays with differential expression, a set of features following a central tendency can be found. In self-self and quasi-self-self experiments, a symmetrical distribution relatively closely following a linear central trend is expected. In step 1808, the method of the present invention seeks to select those features, such as feature 2006, with C1/C2 ratios that follow the characteristic C1/C2 ratio for the data set. The distribution of C1/C2 ratios need not correspond to a linear distribution shown in FIG. 20, but may exhibit a more complex central trend, that might, for example, be mathematically described as a second-order or greater-order polynomial function of the C1/C2 ratios with respect to C 2-signal intensity.

FIGS. 21 A-B and 22-24 illustrate a rank-order method for selecting features, from the set of filtered features produced in step 1806, that follow the central trend in C1/C2 ratios of the features of the data set. In FIG. 21A-B, the C1 and C2-signal intensities for a small set of 25 features is shown. Each feature is described by indices i and j that designate the position of the feature intensity for the feature within each of two 2-dimensional arrays 2102 and 2104 containing feature intensities for a particular channel. Thus, for example, feature (0,0) has the C1-signal intensity 963 (2106 in FIG. 21A) and the C2-signal intensity 1059 (2108 in FIG. 21B).

The C1 and C2-signal intensities for the 25 features can be alternatively considered to reside in two 1-dimensional arrays. FIGS. 22A-B shows the C1 and C2 signal intensities, stored in the 2-dimensional arrays of FIGS. 21 A-B, alternatively considered to be stored in two 1-dimensional arrays, each with a single index f. The feature signals are indexed in row order, with respect to the 2-dimensional arrays, within the 1-dimensional arrays. Thus, the feature signals for feature (1,0) are stored in the sixth elements 2206 and 2209 of the 1-dimensional arrays 2202 and 2204 in FIGS. 22A-B.

The two sets of feature signals within the data set are next ranked with respect to signal intensity. In FIG. 23, two 2-dimensional arrays corresponding to the two 1-dimensional arrays in FIG. 22 contain the C1 and C2 signal intensity data contained in the 1-dimensional arrays in FIG. 22. In the 2-dimensional array 2302, corresponding to 1-dimensional array 2202 in FIG. 22, the first row 2304 contains sorted feature intensities, and the second row 2306 contains the corresponding indices of the feature in the 1-dimensional array 2202 in FIG. 22. For example, the feature with the lowest C1 intensity is contained in the second element 2208 of the 1-dimensional array 2202 in FIG. 22 with index “1,” and is stored the first column 2308 of 2-dimensional array 2302 in FIG. 23, with the entry 2310 in the first row of the first column containing the feature intensity “126” and the entry 2312 in the second row in the first column containing the index “1” of the feature 2208 from the 1-dimensional array 2202 in FIG. 22. The index ρ of the column in a 2-dimensional array (2302 or 2314 ) containing a feature intensity and its associated 1-dimensional array index corresponds to the rank of the feature within the feature-signal data subset stored in the 2-dimensional array and sorted in ascending order with respect to intensity. Thus, feature “1” with intensity “126,” being the lowest intensity feature in the C1-signal data subset, is stored in the first element of array 2302 with ρ=0 and therefore has rank “0.” Array 2302 contains the sorted C1-signal data subset, and the array 2314 contains the sorted C2-signal subset.

In order to find those features having a C1/C2 ratio close to the characteristic C1/C2 ratio for the data set, or, in other words, those features close to the central trend within the data set, the ranks of a feature in the two sorted subsets contained in arrays 2302 and 2314 are compared. Those features having identical or similar ranks in both data subsets are considered to lie within the central trend of the C1/C2 distribution for the data set. For example, feature “1” (2208 in FIG. 22A) has rank “0” in the C1 data subset (2308 in FIG. 23A) and has rank “19” (2316 in FIG. 23B) in the C2 data subset. Feature “1” is therefore not considered to be a central-trend feature. By contrast, feature “21” (2210-2211 in FIGS. 22A-B, 2318 in FIG. 23A, and 2320 in FIG. 23B) has ranks “1” and “0” in the sorted C1 and C2 data subsets, respectively. Feature “21” is therefore considered to be a central-trend feature within the C1/C2 distribution. More specifically, features are considered to be central-trend features when the relative rank displacement is less than a threshold value s, as expressed in the following expression: $\frac{{\rho }_{c1}-{\rho }_{c2}}{\left(\max {\rho }_{c1}-\min {\rho }_{c1}\right)}\le s$ where ρ_C1, is the rank of C1 signal intensity in the sorted C1 data subset, and ρ_C2 is the rank of C2 signal intensity in the sorted C2 data subset. The denominator in the above case can also be considered as the total number of features. In this context, the above formula can be generalized as the normalized relative rank displacement, where the sum of features is the normalization parameter.

The threshold s can be either arbitrarily chosen, or, in an alternative embodiment, is chosen based on the correlation coefficient for the C1/C2 distribution of the features of the data set. The correlation coefficient is provided by the following equation: $\frac{{\left[\sum ^N\left(C1-\stackrel{\_}{C1}\right)\left(C2-\stackrel{\_}{C2}\right)\right]}^2}{\left[\sum ^N{\left(C1-\stackrel{\_}{C1}\right)}^2\right]\left[\sum ^N{\left(C2-\stackrel{\_}{C2}\right)}^2\right]}$ where {overscore (C1)} is the average C1-signal intensity, {overscore (C2)} is the average C2-signal intensity, N is the number of features in the data set. Alternatively, the correlation coefficient can be based on absolute distances from the central trend curve, on vertical distances from the central trend curve, or on any number of other distribution-based density or dispersion metrics.

The selected central-trend features are then sorted according to a combined feature intensity E given by the following expression: E−{square root}{square root over ((C1)²+(C2)²)}

In fact, the features can be sorted by metric of similarity. In case of applying the described mechanism to features with background subtracted intensity, it is necessary to maintain a distinction between features with resultant negative versus positive signals. Therefore, prior to the evaluation of the E metric described above, it is necessary to translate the origin(0,0), such that all features populate the positive quadrant only. In other words, the coordinate axes are moved in order to reposition the ideal, characteristic, low intensity data point, shown in third and fourth quadrant positions in FIG. 17, into the first, positive quadrant. This mechanism is implemented solely to compute E and does not introduce any perturbation into the true data set. FIG. 24 shows the central-trend features selected via the ranking process described with reference to FIGS. 21A-B, 22 A-B, and 23 ranked in ascending order with respect to the combined feature intensity E computed with respect to the above formula. Note that, each feature in the central-trend, combined-feature-intensity data subset stored in the array shown in FIG. 24 is associated with a rank.

In step 1810, the method selects, as chosen features, a lowest combined-feature-intensity portion of the central-trend, combined-feature-intensity features. A threshold-percentile cutoff, x, for choosing the lowest-combined-intensity features may have a specific, fixed value, or maybe tunable with respect to additional constraints, such as the absolute number of features within the chosen percentile range and other constraints. At the end of step 1810, a set of low-combined-intensity, central-trend features have been selected. In one embodiment, the threshold x is determined by starting with an x value that includes nearly all central-trend features, and iteratively reduces x until the slope of the best-fit line and a distribution metric stabilize, assuming that the distribution of data points around the central-trend yields a curve, rather than a line, or alternatively, increases x from a low starting point until the slope of the best-fit line and a distribution metric become non-stable.

In step 1812, a line is fit to the C1/C2 distribution of the low-combined-intensity central-trend features. FIG. 25 illustrates the fitting of a line to the low-combined-intensity, central-trend features. In FIG. 25, the best-fit line 2502 is seen to represent the trend in the distribution of the plotted features, such as 2504. Many different possible line-fitting techniques can be employed. In a preferred embodiment, a technique that minimizes the mean or median absolute deviation (“MAD”) is employed. In one minimization technique, the following quantity M is minimized: $M=\sum _{i=1}^Ny_i-a-b{x}_i$ where y_i is the y coordinate for a plotted feature, x_i is the x coordinate for a plotted feature, b is the slope of the fitted line, and a is the y-intercept for the fitted line.

FIG. 26 illustrates the principle of the above-described minimization technique. FIG. 26 shows the same low-combined-intensity, central-trend features plotted with respect to the fitted line as in FIG. 25. FIG. 26 also shows the vertical distances, such as vertical distance 2602, of plotted features, such as plotted feature 2604, from the best-fit line 2606. The above-described technique minimizes the sum of these vertical distances. Another popular technique for line fitting is the least-squares technique. However, the least-squares technique is not the preferred method, as it is generally more sensitive to outlier data points than MAD-minimization techniques. Other alternatives include median-line fitting techniques. Note that alternative M quantities that can be minimized include the Euclidean distances of data points from the best-fit line, or the median Euclidean distance of data points from the best-fit line.

In step 1814, the line fit to the C1/C2 distribution for the low-combined-intensity, central-trend features is used to further filter this set of features based on their proximity, in the C1/C2 plot, to the best-fit line. FIGS. 27-28 illustrate two methods that can be employed to filter the low-combined-intensity, central-trend features with respect to the fitted line. In a first technique, illustrated in FIG. 27, two threshold lines 2702 and 2704 are constructed parallel with, and above and below, respectively, the fitted line 2706. The distance, in the vertical direction, between the best-fit line 2706 and the threshold lines 2702 and 2704, τ, is related to the aggregate deviation of the plotted features in the C1/C2 plot from the best-fit line, or the value M minimized in the MAD technique described above. The threshold parameter τ may then be calculated by an expression such as: $\tau =\sqrt{\left(\left(\frac{1}{2}\right)*{\left(kM\right)}^2\right)}$ where k is a constant or a tunable parameter. Note that other estimators for the deviation of the plotted features from the fitted line can be substituted for M in the above equation. Note, as well, that other methods for determining the threshold parameter τ may be employed, and that τ itself may be a tunable parameter, rather than calculated from calculated deviations or error metrics for the plotted features. Once the two threshold lines 2702 and 2704 are constructed, it is straightforward to select those low-combined-intensity, central-trend features lying between the two threshold lines 2702 and 2704.

An alternative technique for filtering features based on proximity to the best-fit line is shown in FIG. 28. In this technique, rather than constructing parallel threshold lines, threshold lines with greater 2802 and lesser 2804 slopes than the best-fit line 2806 are constructed to fan out to either side of the best-fit line from the x-axis intercept 2808 or, in other cases, from the y-axis intercept, of the best-fit line. In this case, a threshold angle, τ, 2810 and 2812, can be computed based on a measure of the deviation of plotted features from the fitted line, or can be a tunable parameter. As in the previous threshold technique, those features lying between the threshold lines, or within a cone or hyper-cone, in higher dimensional cases, are chosen. Using either the filtering technique shown in FIG. 27 or that shown in FIG. 28, step 1814 produces a final set of filtered, non-control features.

The final set of filtered, non-control features is augmented, in step 1816, with non-saturated, uniform control features from the full data set, not initially selected in step 1804, that, when plotted in the C1/C2 plot, such as the C1/C2 plots shown in FIGS. 27 and 28, fall between the threshold lines and that also meet the x-percentile cutoff for features ranked with respect to combined intensities. In step 1818, the augmented, final filtered set produced in step 1816 is re-ranked with respect to combined feature intensities, as in step 1808, and a lowest percentile specified by the parameter y, similar to parameter x in step 1810, is selected as a characteristic set of background features. In an optional step 1820, this characteristic set can be further filtered using a box-plot filtering technique in which the highest quartile and/or the lowest quartile features with respect to combined intensities are removed from the characteristic background feature set. As a variation, box-plot analysis with fences is also possible. This essentially widens the scope of the box-plot by including features lying below and above a tunable parameter z times standard deviation of the distribution, in case of a uniform distribution, or z times the interquartile range, with respect to the 2^nd and 3^rd quartiles, respectively. Depending on stringency requirements, this increase in population of the features potentially enhances statistical accuracy.

Next, in step 1822, an ideal background feature data point is constructed from the set of characteristic background features produced in step 1818 or 1820. An ideal background feature is a feature that accurately reflects the offsets in both channels that can be used to background correct a data set. FIGS. 29-30 illustrate the construction of an ideal background data point. In FIG. 29, the characteristic background features are plotted in a C1/C2 plot, distributed about the best-fit line 2902. In FIG. 29, the x,y coordinates for each feature are shown next to the plotted data point, such as the x,y coordinates (9,7) 2904 shown next to the plotted data point 2906. In FIG. 30, an ideal characteristic background feature β is constructed with x,y coordinates equal to the median x coordinate, calculated as the median of the x-coordinate values for the characteristic background features, and the median y coordinate, calculate as the median of the y-coordinate values of the characteristic background features. Thus, in FIG. 30, the ideal background features is plotted 3002 with x,y coordinates (16.5, 7) calculated from the x,y coordinates shown in the table 3004 in the upper right-hand corner of the figure. Note that the coordinates in table 3004 are sorted by x-coordinate value. Sorting by y-coordinate value would more directly reveal the median y-coordinate value 7. Alternatively, the ideal characteristic background feature β coordinates can be calculated as the average x-coordinate and y-coordinate values for the characteristic background features, or may be based on a percentile of characteristic-background-feature x-coordinate values and y-coordinate values or based on the x,y-coordinates of a lowest characteristic background feature.

Next, in step 1824, the ideal characteristic background feature β is further refined, or optimized, by projecting β back to the best fitted line. FIG. 31 illustrates a final refinement technique for calculating the ideal characteristic background feature. In FIG. 31, the data point plotted for the ideal characteristic background feature β 3102 is projected by translating the position of β in a direction perpendicular to the best-fit line 3104 to a position β′ 3106 coincident with the best-fit line 3104. The coordinates for the refined ideal characteristic background feature β′ can be calculated by solving simultaneous linear equations, one for the best-fit line with slope m, and one with slope $\frac{-1}{m}$ that passes through the characteristic background feature β. Although step 1824 may be omitted in alternative embodiments, and the ideal characteristic background feature β used to compute the global, background-signal corrections, the refinement represented by step 1824 has proven to significantly increase the accuracy of the global, background-signal corrections. In optional step 1826, the x coordinate of the ideal characteristic background point is taken as the magnitude of the C2 background-signal correction, and the y coordinate of the ideal characteristic background point is taken as the magnitude of the C1 global background-signal correction, with both the C2 and C1 global background-signal corrections are applied to the entire data set. The net effect for the two-channel background correction is to move an ideal, lowest intensity-ratio data point (1606 in FIG. 16) to the origin of a C1/C2 distribution plot, rather than simply extending a best-fit curve to an x-axis or y-axis intercept, as discussed with reference to FIG. 13. It should be noted that the distance between the unprojected β and projected β′, as well as the difference between the mean and median values of β, may potentially serve as metrics for the quality of feature selection for the 2-channel background correction, and either may serve as the error metric δ returned in step 1824.

In step 1828, the routine “Correct C1 and C2 channel data” is called, in order to accommodate microarray users who attempt to create mock variance stabilization by not subtracting background signal or under-subtracting background signal. Note that failure to subtract the background signal, as described with reference to optional step 1826, may give results that are inconsistent with dye-label-reader bias. Therefore, rather than moving the entire data set, as describe above with reference to optional step 1826, users can move C1 and C2 channel data in accordance with a user-defined constant, referred to as “α.”

FIG. 32 is a control-flow diagram of the routine “Correct C1 and C2 channel data.” The method of correcting C1 and C2 channel data is described with reference to FIG. 32 and concurrent reference to subsequent figures that provide greater detail for many of the steps shown in FIG. 32.

In step 3202, the method of the present invention allows the user to input the user-defined, nonnegative constant α. The user-defined, nonnegative constant α is used to move the C1 and C2 channel data in a manner that is consistent with the dye-label-reader bias. Moving the C1 and C2 channel data in a manner that is consistent with dye-label-reader bias is equivalent to moving the plot of C1/C2 data points parallel to the best-fit line given by: y=b·x+a where b and a are the slope and y-axis intercept of the best-fit line. Typically, the user-defined constant α is employed by the user to ensure that the background correction for C1 and C2 channel data is consistent for inter-microarray as well as intra-microarray hybridization assays.

Next, in step 3204, if the slope b of the best-fit line is greater than “1,” then proceed to step 3206. In step 3206, the value of the user-defined constant α is added to all C2 channel intensity values, and the value α·b is added to all C1 channel intensity values. FIGS. 33A-B illustrate moving a C1/C2 data set having best-fit line slope greater than “1.” In FIG. 33A an entire, hypothetical set of C1/C2 data points, such as data point 3302, having a best-fit line 3301, is plotted. Because the slope b of the best-fit line 3301 is greater than “1,” the entire set of C1/C2 data points are moved with the same direction and magnitude of displacement defined by the vector 3303. The vector 3303 has the same slope b as the best-fit line 3301 with x and y components given by α and α·b, respectively. For example, the data point 3302 is moved parallel to the best-fit line 3301 by adding the x-component a 3304 and the y-component α·b 3305 to the x and y coordinates of the data point 3302, respectively, resulting in a new data point location identified by open circle 3306. FIG. 33B illustrates the net results of moving the entire set of C1/C2 data points by adding the components of vector 3303 to the coordinates of each C1/C2 data point. For example, the characteristic data point 3307 (FIG. 33A) having xy coordinates x′ 3308 and y′ 3309 has been moved to the characteristic data point 3310 (FIG. 33B) having xy coordinates x′+α3311 and y′+α·b 3312.

In step 3204, if the slope b of the best-fit line is less than “1,” then proceed to step 3208. In step 3208, the value of the user-defined constant α is added to all C1 channel intensity values, and the value $\frac{\alpha }{b}$ is added to all C2 channel intensity values. FIGS. 34A-B illustrate moving a C1/C2 data set having a best-fit line slope less than “1.” FIG. 34A is a plot of an entire, hypothetical set of C1/C2 data points, such as data point 3402, having a best-fit line 3401. Because the slope b of the best-fit line 3401 is less than “1,” the entire set of C1/C2 data points are moved by adding the x and y components of vector 3403 to the xy coordinates of each C1/C2 data point. For example, the data point 3402 is moved with the same direction and magnitude as the vector 3403 by adding the x-component $\frac{\alpha }{b}$3404 and the y-component α 3405 to the x and y coordinates, respectively, of the data point 3402 to give the data point identified by the open circle 3406. FIG. 34B illustrates the net results of moving the entire set of C1/C2 data points by adding the components of the vector 3403 to the coordinates of each C1/C2 data point. For example, the characteristic data point 3407 (FIG. 34A) having xy coordinates x′3408 and y′3409 is moved to the characteristic data point 3410 (FIG. 34B) having xy coordinates $x^{\prime }+\frac{\alpha }{b}3411\mathrm{and}{y}^{\prime }+\mathrm{\alpha 3412}.$

In step 3210, if the user is not satisfied with the C1 and C2 data sets after adjustments have been made according to steps 3202-3208, the user has the option of submitting another user-defined constant a and repeating steps 3202-3208. In step 3210, if the C1 and C2 channel data is adjusted to the satisfaction of the user, then the method of the present invention returns to the calling routine “Multi-channel Background Correction.”

The user-defined, nonnegative constant a can be used to adjust C1 and C2 data sets consistently for batches of microarrays, and is applied to the C1 and C2 data set in a dye-bias consistent manner. Note also that, for different feature extraction situations, the method described above with reference to steps 3202-3210 can be applied to background-subtracted as well as non-background subtracted C1 and C2 data.

Adding the offset to both channels as described above with reference to the control flow-diagram in FIG. 32, may significantly improve the dye-normalization process because adjusting the global background to zero alone may leave more features with negative signal in the green channel compared to the red channel (or vice versa). The result is a distortion of the rank consistent data used to perform dye normalization. Adding a user-defined offset, as described above in steps 3206 and 3208, to both channels may bring negative-signal-intensity features back contribute to in the rank consistent normalization method. The results obtained by the method of the present invention are greatly improved if the user-defined offset brings all feature intensities above zero. Ideally, prior to dye normalization, less than 1% of the features should have negative signal intensities after correcting for the background signal according to the method of the present invention.

Although the present invention has been described in terms of a particular embodiment, it is not intended that the invention be limited to this embodiment. Modifications within the spirit of the invention will be apparent to those skilled in the art. For example, the method and system can be straightforwardly extended to determine global background-signal corrections for multi-channel data sets including signals for more than two channels. In one approach, a hyper-dimensional signal-intensity space is considered, each dimension corresponding to a particular channel, and a best-fit hyper-volume calculated for the distribution of feature intensities in the hyper-dimensional signal-intensity space. The lowest-combined-intensity features can then be selected to compute the hyper-dimensional coordinates of an ideal characteristic background data point, from which global background-signal corrections for each channel can be obtained. Alternatively, global background-signal corrections can be calculated iteratively by repeatedly calculating, in a pair-wise fashion, relative global background-signal corrections for pairs of channels. The global background-signal correction method can be implemented in an almost limitless number of ways, in many different computer languages for execution in many different types of computers, using an almost limitless number of different modular organizations, control structures, data structures, and variables. Different sub-methods can be employed. For example, rather than using the rank-ordering method for discovering central-trend features, a geometric or statistical method may be employed. As discussed above, different methods for fitting a line or curve to data-point distributions can be employed. Although each step shown in FIG. 18 is carried out in a preferred embodiment of the present invention, alternative embodiments may omit one or more individual steps shown in FIG. 18, and a workable, background-signal correction method can still obtained, due to the robust nature of the method described in FIG. 18. For example, either or both of steps 1808 and 1824 may be omitted, and the resulting alternative embodiments still produce useful background-signal corrections. Global, background-signal corrections calculated by an embodiment of the method of the present invention can be used to evaluate operation of a microarray reader, calibrate a microarray reader, evaluate the quality of a microarray, evaluate the correctness of background subtraction methods and other data processing methods, and evaluate of the reproducibility of a molecular-array-based experiment. For example, background-signal corrections calculated by an embodiment of the method of the present invention can be compared against standard or expected background-signal corrections, and, if the observed background-signal corrections fall outside the standard or expected ranges, microarrays exceeding or failing to meet standards or expectations may be rejected or remanufactured, microarray readers producing data with observed background-signal corrections falling outside the standard or expected range may be recalibrated, adjusted, or partially remanufactured, and microarray experiments producing data with observed background-signal corrections falling outside the standard or expected range may be redesigned or repeated. The present invention further provides a computer program product that may execute a method of the present invention. The program product includes a computer readable storage medium having a computer program stored thereon and which, when loaded into a programmable processor, provides instructions to the processor of that apparatus such that it will execute the procedures required of it to perform a method of the present invention. The storage medium can, for example, be a portable or fixed computer readable storage medium, whether magnetic, optical or solid state device based.

The foregoing description, for purposes of explanation, used specific nomenclature to provide a thorough understanding of the invention. However, it will be apparent to one skilled in the art that the specific details are not required in order to practice the invention. The foregoing descriptions of specific embodiments of the present invention are presented for purpose of illustration and description. They are not intended to be exhaustive or to limit the invention to the precise forms disclosed. Obviously many modifications and variations are possible in view of the above teachings. The embodiments are shown and described in order to best explain the principles of the invention and its practical applications, to thereby enable others skilled in the art to best utilize the invention and various embodiments with various modifications as are suited to the particular use contemplated. It is intended that the scope of the invention be defined by the following claims and their equivalents:

Claims

1. A method for calculating background corrected signals for a microarray data set having one or more channels, the method comprising: receiving a non-negative constant; selecting a set of low-combined-intensity features from the microarray data set; determining a representation to describe a central-trend of the selected set of low-combined-intensity features within a signal-intensity space; and adjusting the microarray data set based on the non-negative constant and the determined representation.

2. The method of claim 1 wherein selecting the set of low-combined-intensity, features from the data set further includes: selecting non-control features from the data set; filtering the selected non-control features to remove non-uniform features and signal-saturated features; selecting central-trend features from the filtered, selected non-control features and signal-saturated features; and selecting a lowest-intensity percentile subset of the selected set of low-combined-intensity, central-trend features.

3. The method of claim 1 wherein selecting the set of low-combined-intensity, features from the filtered, selected non-control features further includes: ordering each channel-specific data subset within the data set by feature intensity; and selecting as central-trend features those features with identical or similar ranks in all channels.

4. The method of claim 3 wherein selecting a set of low-combined-intensity, central-trend features from the data set further includes: determining a best-fit representation to describe the central-trend of features distributed within the signal-intensity space; and selecting features proximal to the best-fit representation in signal-intensity space.

5. The method of claim 4 wherein selecting a set of low-combined-intensity, central-trend features from the data set further includes: augmenting the selected set of features proximal to the best-fit representation in signal-intensity space with control features of low intensity proximal to the best-fit representation in signal-intensity space.

6. The method of claim 7 wherein determining the representation to describe the central trend of features further includes: constructing a best-fit line, curve, volume, or hyper-volume for two-channel, three-channel, and more-than-three-channel data sets, respectively, that represents the central trend of features distributed within the signal-intensity space.

7. The method of claim 1 further comprising: determining a position of a characteristic background data point based on the selected, low-combined-intensity features in a signal-intensity space with dimensions corresponding to the channels calculating an optional global, background-signal corrections for each channel from the position of a characteristic background data point within the signal-intensity space; and for each channel, selecting the magnitude of the coordinate of the characteristic background data point with respect to the channel in the signal-intensity space as the global, background-signal correction for the channel.

8. The method of claim 7 further includes applying the optional global, background-signal correction for each channel to the data set by subtracting the global, background-signal correction from the feature intensities within the data subset corresponding to the channel.

9. The method of claim 1 wherein adjusting the microarray data set further includes shifting the microarray data set parallel to the determined representation within the signal intensity space.

10. A representation of a background-corrected data set, produced using the method of claim 1, that is maintained for subsequent analysis by one of: storing the representation of the background-corrected data set in a computer-readable medium; and transferring the representation of the background-corrected data set to an intercommunicating entity via electronic signals.

11. Results produced by a molecular-array data processing program employing the method of claim 1 stored in a computer-readable medium.

12. Results produced by a molecular-array data processing program employing the method of claim 1 printed in a human-readable format.

13. Results produced by a molecular-array data processing program employing the method of claim 1 transferred to an intercommunicating entity via electronic signals.

14. A method comprising communicating to a remote location signals which have been background corrected using a global, background-signal intensity correction obtained by a method of claim 1.

15. A method comprising receiving data produced by using the method of claim 1.

16. The method of claim 1 wherein a molecular-array data set having one or more channel includes: a data set containing data subsets corresponding to feature signals obtained from reading a single microarray in two or more different channels; a data set containing data subsets corresponding to feature signals obtained from reading two or more different arrays in a single channel; and a data set containing data subsets corresponding to feature signals obtained from reading two or more different arrays in two or more different channels.

17. Using one or more optional global background-signal corrections calculated by the method of claim 1 to carry out one of: evaluation operation of a microarray reader; evaluation of the quality of background correction; evaluation of the quality of data corrections other than background corrections; calibration a microarray reader; evaluation the quality of a microarray; and evaluation of the reproducibility of a molecular-array-based experiment.

18. A computer program including an implementation of the method of claim 1 stored in a computer readable medium.

19. A method comprising forwarding data produced by using the method of claim 1.

20. A multi-channel, molecular-array data-set processing system comprising: a computer processor; a communications medium by which molecular-array data points are received by the molecular-array-data processing system; one or more memory components that store molecular-array data points; and a program, stored in the one or more memory components and executed by the computer processor, that receives a non-negative constant; selects a set of low-combined-intensity features from the microarray data set; determines a representation to describe a central-trend of the selected set of low-combined-intensity features; and applies the non-negative constant and the representation to correct the microarray data.