"Hook"-calibration of GeneChip-microarrays: Theory and algorithm

The competitive two-species Langmuir model of microarray hybridization

Here the superscript denotes the probe-type (P = PM, MM). The indices "p" and "c" assign probe- and chip-specific parameters, respectively. The probe-index implies the chip specificity as well, i.e. p = p, c. This model predicts that the fraction of "occupied", i.e. dimerized oligonucleotides of a probe spot, Θ_p^P (also called surface probe coverage or occupancy), is directly related to the observed intensity, I_p^P* [14, 18]. The proportionality constant, M_c, specifies the maximum intensity referring to complete occupancy, Θ_p^P = 1, if all oligonucleotides of the respective probe spot on the given chip are dimerized. The minimum intensity referring to the absence of bound transcripts, Θ_p^P = 0, gives rise to the "optical" background intensity, O_c. Throughout the paper we will consider only "net" intensities which have been corrected for the optical background before further analysis, ${\text{I}}_{\text{p}}^{\text{P}}={\text{I}}_{\text{p}}^{\text{P}}\ast -{\text{O}}_{\text{c}}$, using, for example, the zone-algorithm provided by Affymetrix [23].

The latter equation assumes that the large number of different non-specific RNA-fragments in the hybridization solution effectively acts like a single species with the common concentration [N]_c for all probes of the chip [15, 16]. Contrarily, the concentration of specific transcripts, [S]_p, refers to a particular probe sequence, i.e., it represents a "single probe"-property. Microarrays of the GeneChip-type use so-called probe sets of several probes (usually N_set = 11) for estimating the expression of each considered gene. One expects therefore that all probes of a set probe the same, common transcript concentration, i.e. [S]_set = [S]_p for p ∈ set assuming that effects as alternative splicing have been appropriately considered during probe design.

The competitive two-species Langmuir adsorption isotherm (Eq. (1)) considers the effects of non-specific "background" hybridization and of saturation at small and large concentrations of specific transcripts, respectively. The maximum intensity at saturation, M_c, depends on factors such as the number of oligonucleotides per probe spot (which in turn is related to the density of oligomers and to the spot size), the mean number of optical labels per bound target and the settings of the scanner. These factors affect the PM and MM nearly in the same fashion giving rise to virtually identical values of M_c at complete saturation of the probe spots under equilibrium conditions (X_p^P >> 1) [16, 18].

Recent studies report significantly higher limiting intensity values of the PM, compared with that of the MM, i.e. M^PM > M^MM [22]. They interpreted this result assuming a probe-dependent partial dissociation of the duplexes during the post-hybridization washing phase. Another, additional explanation might be the truncation of a considerable amount of the probe oligomers due to incomplete synthesis because this effect causes the asymptote-like flattening of the hybridization isotherms at intermediate and large transcript concentrations in a sequence-dependent manner [1, 9].

We will apply in the following analysis the special-case of the common intensity asymptote for all probes of the chip according to Eq. (1). Possible consequences of deviations from this assumption for the data analysis will be addressed in a separate study.

Matched and mismatched microarray probes

This S/N-ratio, R, provides the specific binding strength of the PM in units of the non-specific one. It can serve as a relative measure of the expression degree because it is directly related to the concentration of specific transcripts, [S]_p. It scales the expression degree in a probe-specific fashion.

Part a of Figure 1 shows the courses of the intensities of a typical PM/MM pair as a function of the parameter R (see Eq. (6)). The PM intensity sigmoidally increases from its minimum value, I_p(R = 0), to I_p(R = ∞) = M_c, at small and large abscissa values, respectively. The respective probes referring to these limiting cases are either exclusively non-specifically hybridized or completely saturated with surface coverages of Θ_p^PM(0) = X_p^PM,N/(1 + X_p^PM,N) ≈ X_p^PM,N and Θ_p^PM(∞) = 1, respectively. The concentration and S/N-ratio referring to the inflection point of the isotherm at half-way between these values are

respectively. They specify the condition at which 50% of the free probes available in the absence of specific transcripts become occupied. The approximations at the right-hand side of Eq. (7) refer to small X_p^PM,N << 1.

which shows that the horizontal shift between the PM- and MM-isotherms is log(s_p) and log(s_p/n_p) in the log-scale of [S] and R, respectively.

The delta- and sigma-transformations

The MM probes were designed as reference for estimating the non-specific background contribution to the respective PM intensity [2, 5, 25]. The "simple" subtraction of the MM-intensity from that of the PM however partly failed as correction method because both probes differently respond to non-specific and specific hybridization due to their complementary middle bases which, for example, gives rise to negative PM-MM intensity differences [15].

According to the Langmuir model, the behaviour of the PM and MM can be understood on the basis of the same hybridization rules where both probe types however differ with respect to their effective association constants for probe/target dimerization (see above). The intensities of the PM and MM are consequently expected to correlate in a well defined fashion. This mutual relation is determined by the mismatch design of the reference probe, the particular probe sequences and by the concentrations of specific and of non-specific RNA fragments in the sample solution used for hybridization on the particular chip [18].

In summary, the hyperbolic intensity functions of the PM and MM can be transformed into Δ and Σ coordinates which are governed by essentially four parameters, the start values Δ_p^start ≅ Δ_p(0) and Σ_p^start ≅ Σ_p(0) and the exponents α_p and β_p. They were chosen to provide a simple geometrical interpretation of the Δ-vs-Σ trajectory in terms of its start-coordinates and its vertical and horizontal dimension with respect to the start values (see below and Figure 1 and Figure 2).

The hybridization regimes

Part b and c of Figure 1 show the transformed intensities taken from part a of the figure as a function of the parameter R = R_p^PM. The course of the log-difference, Δ_p(R), can be roughly divided into five regimes which reflect different hybridization characteristics of the PM and the MM probes with increasing degree of specific hybridization (see Figure 1, part b):

(1) N-regime: In the non-specific-regime, at small values R → 0, both, the PM and MM nearly exclusively hybridize with non-specific transcripts. Saturation can be typically neglected in this range (B_p^P ≈ 1, see Eqs. (10) and (11)). The limiting ordinate value for X_p^P,N << 1 estimates the ratio of the binding constants referring to the respective pair of complementary middle bases in the PM and MM sequences, Δ_p(0) ≈ log n_p (see Eq. (4)). We will use the approximation of weak non-specific binding throughout the paper.

(2) mix-regime: In the subsequent mixed-regime, both, specific and non-specific transcripts significantly contribute to the observed intensity of the probes. The log-difference Δ increases with increasing amount of specific transcripts. The positive slope of Δ_p(R) implies ∂Δ/∂R ~ (1-10^-α) > 0, and thus α_p > 0 or equivalently s_p > n_p (see Eq. (12)). The increase of Δ_p(R) consequently reflects the simple fact that the specific binding constant of the PM exceeds that of the respective MM, i.e., K_p^PM,S > K_p^MM,S, if one assumes K_p^PM,N ≈ K_p^MM,N (see below and Eq. (4)).

This rough approximation assumes Δ_p(0) << 1 <β_p and R_max >> 1. At conditions of weak saturation Eq. (13) simplifies with 0.5 (β_p - α_p) >> 1 into ${\Delta }_p\left(R_{\max }\right)\approx {\Delta }_p^{Linear}\left(\infty \right)={\Delta }_p(0)+{\alpha }_p>0$. At these conditions the height of the maximum directly provides the log-transformed PM/MM-ratio of the specific binding constants, α_p.

(4) sat-regime: In the saturation-regime the probes become progressively saturated with bound transcripts (B_p^P > 1). This effect first and foremost affects the PM due to their higher specific binding constant (see above). As a consequence Δ_p starts to decrease.

(5) as-regime: At very large expression degrees both, PM and MM reach their maximum intensity upon complete saturation. In this asymptotic-regime the trajectory reaches the abscissa for R → ∞, Δ_p(∞) ≈ 0 (see Eq. (10)).

The respective log-sum of the intensities, Σ_p(R), is shown in part c of Figure 1. It varies in a similar, sigmoidal fashion as the individual log-intensities of the PM and MM (compare part a and c of Figure 1). Here the mix-, S- and sat-regimes merge into one region of increasing Σ whereas the N- and as-regimes provide the minimum and maximum values, ${\Sigma }_{\text{p}}(0)\approx {\Sigma }_{\text{p}}^{\text{start}}$ and Σ_p(∞) = log M_c, respectively. With Eqs. (11) and (12) one obtains for the difference

β_p ≈ Σ_p(∞) - Σ_p(0) > 0

Eqs. (15) and (14) provide ${\Sigma }_{\text{p}}\left({\text{R}}_{\max }\right)-{\Sigma }_{\text{p}}(0)\approx \frac{1}{2}\left({\Sigma }_{\text{p}}\left(\infty \right)-{\Sigma }_{\text{p}}(0)\right)>0$, i.e., the maximum of Δ _p (R) roughly bisects the total range of the Σ-coordinate.

The Δ-vs-Σ trajectory

In the next step we plot the transformed intensities into Δ-vs-Σ coordinates (see Figure 2, part a). This presentation, also known as M-vs-A plot (difference-vs-sum), reflects the binding isotherms of a PM/MM-probe pair. The obtained Δ-vs-Σ trajectory shows a characteristic curved shape with start-, end- and maximum-points referring to the S/N-ratios R = 0, R = ∞ and R = R_max, respectively. They consequently define the N-, as- and S-hybridization regimes. The mix- and sat-regimes can be attributed to the increasing and decreasing parts of the Δ-vs-Σ trajectory, respectively.

The parameters α_p and β_p define the height and the width of the obtained Δ-vs-Σ curve (see also Eqs. (13) and (14)). The Δ- and Σ-coordinates of the characteristic points depend on the PM/MM-ratios of the binding constants (see Eq. (4)), on the maximum intensity, Σ_p(∞) = logM_c, and on the mean intensity of the chemical background due to non-specific hybridization, Σ_p(0) ∝ log(I_p^PM(0)) + log(I_p^MM(0)). Hence, the Δ-vs-Σ trajectory links the observed probe intensities with essential hybridization characteristics in terms of simple geometric parameters.

The horizontal scale of the Δ-vs-Σ trajectory

Eqs. (16) and (17) provide mean values averaged over the respective PM/MM-probe pair. The "individual" coverages of the PM and MM probes, Θ_p^PM and Θ_p^MM, and the respective fraction of specifically hybridized oligomers, x_p^PM,S and x_p^MM,S, in addition depend on the relative Δ-coordinates Δ-Δ(∞) and Δ-Δ(0), respectively (see Eqs. (42) and (45) in the Appendix A).

Part b of Figure 2 shows the surface coverage and the fraction of specifically occupied oligomers for the Δ-vs-Σ trajectory plotted in part a of the figure. Note that x^{P, S} and Θ^P exponentially scale with the coordinate differences Σ-Σ(0) and Σ-Σ(∞), respectively (see Eqs. (17) and (16), respectively).

Consequently, the fraction of specifically occupied probes steeply increases in the raising part of the Δ-vs-Σ trajectory (mix-regime) whereas the probe coverage steeply increases in its decaying part (sat-regime). The contribution of non-specific hybridization and/or the effect of saturation of a particular probe can be essentially neglected if the distance of its Σ-coordinate from the start and/or end points exceeds unity. Particularly, one obtains <x_p^S> > 0.9 for Σ-Σ(0) > 1 and < Θ_p> < 0.1 for Σ(∞)-Σ < 1.

The horizontal shift between the PM-and MM-curves in part b of Figure 2 illustrates the "delayed response" of the MM with respect to the specific transcript concentration: The MM reach a certain ordinate-level of the surface coverage and of the fraction of specifically bound probes at larger abscissa values and thus at larger concentrations of specific transcript concentrations than the PM (see also Eqs. (7) and (8)).

For abscissa values Σ < Σ(∞) -1, Eq. (18) simplifies into log(<R> + 1) ≈ Σ-Σ(0). Hence, the Σ-axis nearly linearly scales with the logarithm of the mean S/N-ratio. For the S/N-ratio of the PM, this equation modifies into $\log ({\text{R}}_{\text{p}}^{\text{PM}}+1)\approx \left(\Sigma -\Sigma (0)\right)+\frac{1}{2}\left(\Delta -\Delta (0)\right)$ (see Eq. (46) below), i.e., it depends in addition on the vertical coordinate of the Δ-vs-Σ trajectory.

In summary, the position of a probe-point along the Σ-coordinate estimates the hybridization degree of the respective probe spot in terms of relative concentration measures characterizing either the S/N-ratio (Eq. (18)), the relative occupancy of the probe oligomers with specific transcripts (Eq. (17)), their overall degree of occupancy of (Eq. (16)) and the specific binding strength of the considered probe pair (Eq. (19)).

The probe coverage (Eq. (16)) provides an additional interpretation of the horizontal dimensions of the Δ-vs-Σ trajectory: For the N-point one obtains with Σ = Σ(0) the coverage due to non-specific hybridization, $〈{\Theta }_{\text{p}}^{\text{N}}(\text{R}=0)〉\approx {10}^{\left({\Sigma }_{\text{p}}(0)-{\Sigma }_{\text{p}}(\infty )\right)}={10}^{-{\beta }_{\text{p}}}$, because (almost) exclusively non-specific transcripts bind to the probes. Note that this "non-specific" coverage is exponentially related to the "width"-exponent, β_p (see Eq. (14)), and thus to the horizontal distance between the N- and the as-points. The remaining, not-occupied and thus free oligomers serve as potential binding sites for specific targets, i.e., $〈{\Theta }_{\text{p}}^{\text{free}}(\text{R}=0)〉=1-{10}^{-{\beta }_{\text{p}}}$. The horizontal dimension of the Δ-vs-Σ trajectory consequently specifies the maximum amount of free probes available for specific binding at R = 0 and thus the measurement range of the probe spots for estimating the expression degree. The narrowing of the model curves reflects the diminishing capacity of the respective probes for specific transcript binding. Figure 3 (part a) illustrates the narrowing of the Δ-vs-Σ trajectory upon increasing the non-specific background contribution. The special ideal case β = -∞ consequently refers to hybridization without non-specific background.

The vertical scale of the Δ-vs-Σ trajectory

The discrimination score roughly estimates the fraction of the signal due to specific hybridization (see [26]). The approximation on the right hand side of Eq. (20) seems save for small values DS << 1.

The discrimination score serves as the basic parameter in the MAS5-algorithm to calculate the so-called detection call (DC) which judges the "presence" or "absence" of a gene. Hence, the vertical scale of the Δ-vs-Σ trajectory is related to the detection call: the higher the Δ_p-value of a probe the higher the probability of the presence of the respective specific transcript in the hybridization solution. We will discuss this point more in detail in the accompanying paper in connection with our alternative method for classifying the genes into present and absent ones (see below).

Here Δε₁₃^WC-SC denotes the dimensionless free energy gain (given in units of the thermal energy, RT) upon replacements of the type B•b^c → B^c•b^c (i.e. WC → SC) for the base B_p = A, T, G, C at sequence position 13 of the probe (for example, C•g → G•g; upper case letters refer to the DNA-probe; lower case letters refer to the bound RNA-fragment, b = a, u, g, c; the superscript "c" indicates the respective complement). Accordingly, Δε₁₃^WC-WC is the respective free energy change upon WC-reversals, B•b^c → B^c•b (for example, C•g → G•c).

Hence, the ordinate position of the starting point of the Δ-vs-Σ trajectory estimates the effective free energy change upon replacing the central base in complementary WC-pairings, i.e. Δ_p(0) ≈ -Δε^WC-WC(B_p) (see Eqs. (11) and (21)). The relative ordinate value of the maximum is related to the respective free energy change upon replacing the central WC-pairing in the specific PM-duplexes by the respective SC-pairing in the MM-duplexes, i.e. Δ_p(R_max) ≈ -Δε^WC-SC.

Figure 3 illustrates that the maximum height of the Δ-vs-Σ trajectory starts to decrease for relatively small widths referring to large strengths of non-specific hybridization (β_p < 3) because saturation onsets almost in the mix-range. In such cases the observed vertical dimension of the trajectory potentially underestimates the height-parameter α_p (see Eq. (13)) which however can be obtained by appropriate curve fitting using Eq. (10) (see below).

In summary, the Δ-vs-Σ trajectory spans a sort of natural or intrinsic metric system between distinctive points which characterizes the binding thermodynamics of the probes of the particular microarray. The horizontal dimension characterizes the measurement range of the respective probe whereas the vertical dimension reflects the free energy gain due to the change of the central base pairing in the respective duplexes of the PM and MM probes.

Δ-vs-Σ trajectories of individual probes

Each probe is characterized by its "individual" Δ-vs-Σ trajectory which describes the intensity change upon increasing content of S-transcripts in the range 0 ≤ R ≤ ∞. We used the Affymetrix HG-133 spiked-in data-set to study the R-dependence of selected probes http://www.affymetrix.com/support/technical/sample_data/datasets.affx. This data set was generated by Affymetrix to calibrate the observed intensities on the basis of known transcript concentrations. Particularly, transcripts referring to 42 selected genes were titrated with increased concentration onto a series of chips using the Latin-squares design. The non-specific background was taken into account by adding a HeLa-cell line extract to all hybridizations which does not contain the spiked-in transcripts.

Part a of Figure 4 shows the trajectories of six selected probes together with fits by means of Eq. (10) (compare curves and symbols). The probe-labels 1 to 6 are chosen to increase with increasing number of C and decreasing number of A per probe sequence (see Figure 4). The observed intensities and thus also the trajectories are functions of the binding constants for DNA/RNA duplex formation, which in turn depend on the sequences of the 25 meric probes. For example, the binding affinity of C•g WC-pairings exceeds that of A•u pairs in the hybrid duplexes. In general, the probes with a higher amount of cytosines are therefore expected to bind the RNA fragments more strongly than probes with a higher amount of adenines. Equation (3) predicts for the increase of K_p^P,N (and of logX_p^PM,N, see Eq. (3)) the decrease of the horizontal dimensions, β_p, of the respective probe-trajectory.

Indeed, the increase of the cytosine-content causes the narrowing of the trajectories by the shifting of their start-point, Σ_p(0), towards larger abscissa values at invariant Σ_p(∞) = const., which is assumed to be constant across all probes because of their common maximum binding capacity. Note that the width of the trajectories and thus the binding strength of the non-specific background varies over about two orders of magnitude, logX_p^PM,N ≈ -4 to -2, for the six selected probes.

The Δ-coordinates of the starting- and maximum-points of the selected probe-trajectories show considerable variation without obvious correlation to their sequence characteristics. We calculated the trajectories of all spiked-in probes (~500) using the results of our previous analysis of the hybridization isotherms (see refs. [18] and [16] for details) to estimate the variance of the positions of their starting- and maximum-points. The boxplot in part b of Figure 4 visualizes the center and the width of the distributions of the obtained Δ_p(0)- and Δ_p(R_max)-data in vertical and horizontal directions.

The respective coordinates of the individual probe-trajectories depend mainly on the particular probe pairings of the middle bases in the non-specific and specific duplexes, respectively (see Eqs. (11) – (13) and (21)). To filter out the underlying sequence effects we calculated "mean" trajectories for all probe pairs with a certain middle base (see Figure 4, part b). These middle-base related mean trajectories are shifted each to another in vertical direction according to C ≈ T > G ≈ A for the N-, and C > G ≈ T > A for the S-point, respectively. This systematic trend is in agreement with Eq. (21) which predicts that the vertical positions of the N- and S-points are functions of the middle base of the respective probe sequences. The observed relations reflect the purine-pyrimidine asymmetry of binding strength of complementary WC-pairings at the N-point (i.e., Δε₁₃^WC-WC(B) ≠ Δε₁₃^WC-WC(B^c)) and the higher stability of the WC-pairings compared with SC-mismatches at the S-point, Δε₁₃^WC-SC(B) (see Eq. (21)). Note that the specific binding constants of PM exceed that of the MM on the average by the factor of s ≈ 7 whereas for non-specific binding one obtains a mean ratio of n ≈ 1.2.

The comparison of the middle-base related trajectories with the width of the N- and S-boxes indicates that the systematic effect due to the middle-base explains the variability of the Δ_p(0)- and Δ_p(R_max)-coordinates in the limits of their 25% and 75% quartiles. The consideration of the nearest neighbors of the middle base further broadens this range: For illustration we show the respective mean trajectories for the "middle triples" CCC and CGC which provide the strongest and weakest binding affinities among the 64 possible combinations of three adjacent bases, respectively (see [18] and [13]).

In summary, the transformed intensity data of individual probes are well described by the Δ-vs-Σ trajectories predicted by the Langmuir-isotherms. The presented data illustrate the probe-specific variability of the Δ-vs-Σ trajectories due to sequence effects. The positions of the start- and maximum-points can be attributed to the differences between the PM and MM probe-sequences which affect the respective binding constants in a middle-base dependent fashion.

The raw hook curve

To construct the mean Δ-vs-Σ trajectory we plot the PM-MM log-intensity difference of all probe pairs of a particular chip, Δ_p, versus the set-averaged log intensity of the respective probe set, <Σ>_set, in a first step (see Eq. (9) and part a of Figure 5). Additional set-averaging of the log-difference, <Δ>_set, reduces the scatter width of the data points along the ordinate roughly by the factor of ~ √N_set ~ 3 (see the yellow dots in Figure 5). In the next step we smoothed these data by calculating the moving average over a sliding window of N_mov ≈ 100 subsequent probe sets along the abscissa to extract the mutual dependence between <Δ>_set and <Σ>_set. The resulting plot is called (raw-) hook-curve because of its typical shape (see part b of Figure 5). Each probe-set is characterized by its "hook"-coordinates, Σ^hook = <Σ>_p∈set and Δ^hook = <<Δ_p>_p∈set>_mov.

The shape of the hook curve basically agrees with that of the Δ-vs-Σ trajectory (compare with Figure 2). We attribute the rising and decaying part and the maximum in-between to the mix-, sat- and S-regimes, respectively, which refer to the mean hybridization level of the respective probes on the chip. The hook-curve obviously does not reach the asymptotic as-regime with Δ(∞) ≈ 0. This result does not surprise because complete saturation is usually circumvented by reasonably chosen hybridization conditions. Otherwise the probes completely lose their sensitivity to detect changes of the transcript concentration [14].

At small abscissa values the hook curve starts with a virtually horizontal part (slope(N) < 0.1) which is separated from the subsequent mix-range (slope(mix) > 0.6) by a distinct break-point. In the appendix we show that the observed initial tiny slope can be explained by the relative strong correlation between the intensities of non-specifically hybridized PM- and MM-probes in the absence of specific transcripts (ρ > 0.7, see Eq. (49)) which in turn seems reasonable in view of the nearly equal sequences of both probe-types which suggest similar variations of the sequence-specific bindings strengths from probe pair to probe pair. On the other hand, for the mix-region the hybridization model predicts a considerably increased slope of slope(mix) ≈ 1.15·α_p ≈ 0.9 > slope(N) ≈ 0.8 (see Eqs. (47) and (49) with·α > 0.8 and ρ > 0.7, respectively). Hence, the break in the course of the hook-curve can be explained by the onset of specific hybridization for probes (and probe-sets) located on the right from this point.

Accordingly, the break-point between the N- and mix-ranges was used to classify the probe-sets into two sub-ensembles: (i) "absent"-ones for Σ^hook < Σ^break, which are assumed to hybridize predominantly non-specifically owing to the absence of specific transcripts in the hybridization solution (R = 0); and (ii) "present"-ones for Σ^hook > Σ^break, which significantly hybridize specifically owing to the presence of specific transcripts (R > 0).

The exact position of the break point of a particular hook curve was estimated by a simple algorithm based on the joint least-squared fit of the Δ^hook-data to a linear function for Σ^hook < Σ^break, and a quadratic function for Σ^hook > Σ^break (see Figure 5). The algorithm systematically varies the position of the break, Σ^break finally returning the optimum value of Σ^break which best fits the data.

Sensitivity-corrected intensity-data and sensitivity profiles

where I^P_0,p denotes the corrected intensity of probe p. The sequence-specific incremental contribution, δ A_p^P (P = PM, MM), the so-called sensitivity of the probe, additively splits into two terms due to specific and non-specific hybridization, δ A_p^P,N and δ A_p^P,S, which are weighted by the fraction of non-specifically and specifically hybridized oligomers, $x_p^{P,N}=\min (L_p^P(0)/L_p^P,1)$ and x_p^P,S = (1 - x_p^P,N) with logL_p^P(0) ≈ logI_p^P(0) = <logI_p^P>_p∈N + δ A_p^P,N and L_p^P = I_p^P/(1-I_p^P/M_c), respectively. The brackets <...>_p∈N denote averaging over all probes from the N-range of the hook curve.

with δ_k(b_m, ξ^P_k,m) = 1 for b_m = ξ^P_k,m and δ_k(b_m, ξ^P_k,m) = 0, otherwise. Here ξ^P_k,m denotes a subsequence of m consecutive bases starting at position k of the respective probe sequence. Each set of sensitivity terms consequently comprises 25 × 4 = 100, 24 × 16 = 384 and 23 × 64 = 1472 parameter values for the SN, NN and NNN models, respectively.

Altogether, the intensity-correction requires four sets of such profiles, δε_k^P,h(b_m), referring to P = PM, MM and h = N, S, respectively. We used weighted multiple linear regression of the normalized intensity data taken from selected sub-ensembles of probe sets to estimate the required sensitivity-profiles (see Appendix C and [18]). These sub-ensembles were chosen to comprise probes which are predominantly hybridized with non-specific (p ∈ N) or specific (p ∈ S) transcripts for estimating δε_k^P,N(b_m) and δε_k^P,S(b_m), respectively. Probes of the former sub-ensemble are taken from the initial horizontal part of the hook curve which was assigned to the N-regime meeting the condition Σ_p^hook < Σ^break (see Figure 5). The respective absent-probes are assumed to bind predominantly non-specific transcripts.

The second sub-ensemble of predominantly specifically hybridized probes (p ∈ S) was selected according to the condition Σ_set^hook > Σ^80% where Σ^80% defines a threshold referring to <x^S> > 0.8, i.e. to probe spots with a fraction of specifically-hybridized oligomers of at minimum 80%. According to Eq. (17) the respective probes are selected to meet the "horizontal distances"-criterion with respect to the break-point (Σ^80% - Σ^break) > -log(0.2) ≈ 0.7. The somewhat arbitrary choice of <x^S> > 0.8 turned out to be not very crucial with respect to the quality of the obtained correction. On one hand the higher the value of <x^S> the smaller the residual contribution of non-specific hybridization in the selected sub-ensemble and the better the obtained sensitivity profiles characterize specific binding [16]. On the other hand, the number of probe sets in the sub-ensemble decreases with increasing <x^S> giving rise to more noisy sensitivity profiles and thus to less precise correction terms [16]. The chosen value of <x^S> >0.8 provides a good compromise (see Figure 5).

Typical SN- and NN-sensitivity profiles of the PM and MM referring to the N- and S-subsets are shown in Figure 6. Note their different shape, especially that of the S-profiles of the MM with the typical "dent" in the middle of the sequence. It is caused by the mismatched base pairing in the specific duplexes of the MM which give rise to changed interaction characteristics compared with WC-base pairings. We refer to our previous papers for the detailed discussion of the sensitivity profiles in terms of base-pair interactions in the respective probe-target duplexes [15, 16, 18]. In the present context it is important to discriminate between the four different sensitivity profiles to properly correct the intensity data for sequence specific effects.

Importantly, for perfect selection of the N-subset one expects the same sensitivity profiles for the PM and MM probes, since the nonspecific background bears no particular relationship to any of both kinds of probes. The profiles shown in Figure 6 confirm this prediction. The degree of similarity of the N-profiles of the PM and MM probes thus provides a criterion for the proper selection of the N-subsets.

The corrected hook curve

In the next step, the corrected intensity values were transformed into (Δ_0,p, Σ_0,p)-coordinates using the corrected intensities in Eq. (9). Subsequently the (Δ_0,p, Σ_0,p)-data were set-averaged and smoothed in the same fashion as described above for the non-corrected intensities to get the sensitivity-corrected version of the hook-curve with the coordinates (Σ₀^hook, Δ₀^hook).

The sensitivity-correction reduces the scattering of the probe- and probe-set-level intensity data about the smoothed hook-curve (see Figure 7). With the N-, mix-, S- and sat-hybridization regimes it essentially shows the same features as the uncorrected hook curve. The break between the N- and mix-regimes was again used to classify the probe sets into absent and present ones. The sensitivity correction affects the (Σ₀^hook, Δ₀^hook)-coordinates of each probe set with possible consequences for its classification into absent and present ones. In a second iteration we therefore re-calculated the four sensitivity profiles on the basis of corrected hook-curve by applying the same criteria for probe selection as above, i.e., Σ₀^hook < Σ₀^break and Σ₀^hook > Σ₀^80% for the N- and S-profiles, respectively. Typically the re-calculated sensitivity profiles only marginally change compared with the profiles which were obtained on the basis of the uncorrected hook curve indicating the relative robustness of the chosen classification criteria.

The detailed comparison between both versions of the hook-curve reveals subtle differences especially of the initial N-regime: After correction it changes slope and, more importantly, considerably narrows in horizontal direction by roughly 50% (compare Figure 7 and Figure 5). The steeper slope of the N-region after correction reflects the reduced cross correlation between the corrected PM- and MM-probe level data (see appendix B). Typically, the coefficient of correlation decreases from values ρ ≈ 0.7–0.9 of the uncorrected data to ρ ≈ 0.4–0.7 after correction which gives rise to the change of the slope by the factor of 1.5 – 3 (see Eq. (49)).

The narrowing of the N-range results from the partial removal of sequence effects which, to a large degree, cause the variability of probe signals [14]. Ideally, the correction of the intensities for sequence specific effects is expected to shrink the N-region into one point referring to one and the same corrected background intensity for all probes of the chip (see Eqs. (11) and (3)). The residual width of the N-region reflects the deficiency of the method due to at least three, essentially undistinguishable effects: Firstly, the "fit"-error of the position-dependent sensitivity model which only incompletely corrects the intensity data for sequence effects caused, e.g. by the specific folding of probe or target and/or bulk dimerization between different targets [1]; secondly, the "N-concentration" error due the simple assumption to consider all non-specific transcripts as one effective species (Eq. (3)) and thirdly, the "classification" error of the simple break-criterion which imperfectly distinguishes between "present" and "absent" probes due to its limited specificity.

Fit of the hybridization model

The corrected hook curve manifests the mean hybridization characteristics of the probes on the particular chip in sensitivity-corrected (Σ₀^hook, Δ₀^hook)-coordinates. It shows essentially the same properties as the Σ-vs-Δ trajectory of a single probe except the N-regime. For quantitative analysis we fit the hybridization model introduced above (Eq. (10)) to the corrected hook data in the mix-, S-, sat- and as-ranges (i.e. at Σ₀^hook > Σ₀^break, see Figure 8). The least-squares gradient descent algorithm searches for optimum values of α_c, β_c and Σ_c(0) ≈ Σ₀^break and Δ_c(0) ≈ <Σ₀^hook>_p∈N. Note that the break position Σ₀^break slightly deviates from the centre of the N-range <Σ₀^hook>_p∈N because of the width of the N-regime (see Figure 7). We define the total width of the hook curve as ${\beta }_c={\Sigma }_c\left(\infty \right)-{〈{\Sigma }_0^{hook}〉}_{p\in N}$.

scales the concentration ratio of specific and non-specific transcripts with the respective ratio of chip-averaged mean binding constants of the PM for specific and non-specific binding. Note that this scaling is common for all probes of the chip in contrast to the probe-specific scaling of R_p^PM in Eq. (6).

In summary, the hook-curve analysis of the intensity data thus provides chip characteristics such as the mean contribution of the non-specific background, the asymptotic intensity maximum and the mean PM/MM-sensitivity gain, and, in addition, the expression degree on probe-set basis in terms of the S/N-ratio.

Signal distributions

Part a of Figure 9 shows the probability-density distribution to find a probe pair a at a given position of the hook-abscissa, p(Σ₀^hook) = ΔN/(N_tot·ΔΣ^hook·) (here ΔN/N_tot is the fraction of probe-pairs found per abscissa interval ΔΣ^hook). About 40% of the probes fall in the N-range in this example and thus are classified "absent". We separately plot the density distribution of these absent probes as a function of their corrected log-intensities (part b of Figure 9). The obtained PM- and MM-probe-level data are well described by normal distributions (P_p^N(x, σ_c) = N(μ_c, σ_c) with $x\equiv x_p={\left(\log I_{0,p}^P-{\mu }_c^P\right)}_{p\in N}$) of with virtually the same width (σ^PM ≈ σ^MM) and slightly shifted centres, ${\mu }_{\text{c}}^{\text{PM}}>{\mu }_{\text{c}}^{\text{MM}}({\mu }_c^P={〈\log I_{0,p}^P〉}_{p\in N})$.

The total distribution decays exponentially to a good approximation at Σ₀^hook > Σ_c(0) in many cases (part a of Figure 9). Define the parameter $r\equiv \log X_{set}^{PM}-{\mu }^{PM}=\log (R+1)$, related to the S/N ratio R.

The exponential behaviour is approximately preserved with respect to this parameter in the mixed and S ranges since $r\approx \left({\Sigma }^{hook}-{\Sigma }_c(0)\right)+\frac{1}{2}\left({\Delta }^{hook}-{\Delta }_c(0)\right)$ from Eq. (18) (see part c of Figure 9). Note that r directly relates to the binding strength of specific hybridization in the relevant range of intermediate and large R-values, r ~ logR ~ log(X_set^PM,S). The initial value of the distribution at r = R = 0 provides the fraction of virtually absent probes, f^absent, whereas the area under the normalized distribution for r > 0 consequently gives the fraction of present probes, f^present = (1 - f^absent).

In the intermediate range the signal is approximated by a constant with f₂^present = 1 - (f^absent + f₁^present) (f^present = f₁^present + f₂^present) whereas for r = 0 one gets the fraction of absent probes. The decay constant λ_r defines the r-range which decreases the probability of the specific signal by one order of magnitude. It thus defines a characteristic S/N-ratio (in logarithmic scale) of the chip which characterizes log-ratio of the mean specific and non-specific signals, or in other words, to which extend the specific signal exceeds the non-specific one. The S(pecific)/N(onspecific)-ratio thus can be also interpreted as a sort of S(ignal)/N(oise)-ratio of a given probe set. The value of the r-related decay constant λ_r slightly exceeds that of the hook-related distribution owing to the larger range covered by the PMonly S/N-ratio (λ_r = 1.0 versus λ = λ_Σ = 0.75; compare part c and a of Figure 9).

An estimate of the mean decay constant can be obtained by simple averaging over the respective R-range

⟨λ⟩ ≈ (⟨log(R + 1)⟩_{R > R'}·In 10)

Note that the exponential decay in Eq. (28) is equivalent with the power law, 10^{-r/λ r}= (R+1)^{-1/λ r}, which has been shown to describe the probability distribution of expression values in a series of microarray experiments [28–31]. This power-law function is known as the Zipf's law, observed in many natural and social phenomena. The presence of such power-law function in principle prevents an intrinsic cut off point between "on" genes and "off" genes. The analysis of expression data in terms of their probability distribution therefore is of basic importance for judging the expression level in terms of global characteristics. The hook-transformation provides such characteristics in terms of the signal distribution along the Σ^hook and/or r coordinates for the mix-, S- and sat-hybridization ranges.

Let us now consider the S/N-ratio of a particular probe at two conditions: (i) with x = 0, referring to the non-specific binding strength in the centre of the normal background distribution, X_p^PM,N ≈ 10 ^μ /M_c:$R_p=M_c\cdot X_p^{PM,S}/{10}^{\mu }$; and (ii) with x ≠ 0, i.e. X_p^PM,N(x) ≈ 10^μ+x/M_c and $R(x)=M_c\cdot X_p^{PM,S}/{10}^{x+\mu }$.

After transformation into logarithmic scale and combination of both situations we get log R(x) = log R _p - x and for R > 1 to a good approximation r(x) ≈ r - x.

Note that the upper integration range is effectively restricted to x ≤ r because of p^S(r-x,λ) = 0 for x > r (see Eq. (28)). Eqs. (30)–(31) define the convolution product of the background and specific-signal distributions. A similar approach was used in ref. [30] and partly also in the RMA-algorithm [29].

De-saturation of the signal and convolution based expression estimates

Figure 10 re-plots the hook curve shown in Figure 8 after linearization of the corrected intensities using Eq. (32). The sat-range now levels off into the asymptotic value α_c (compare with Figure 8). This "de-saturated" signal additively decomposes into contributions due to non-specific and specific hybridization (Eq. (2)). We now calculate the expression value as the weighted "glog"-mean of the specific signal,

Particularly, Eq. (33) corrects the linearized signal (Eq. (32)) for the probe-specific background contribution using the mean background of the chip and the probe-specific sensitivity for non-specific hybridization, $L_p^{P,N}(x,{\mu }_c^P)={10}^{\delta A_p^{P,N}+{\mu }_c^P+x}$ (Eq. (23)). The use of the generalized logarithm, $\text{g}\log (\text{x})\equiv \log \left(\frac{1}{2}(\text{x}+\sqrt{{\text{x}}^2+\text{c}})\right)$, for data transformation is advantageous in two respects: Firstly, it enables processing negative arguments which may appear due to data scattering and, secondly, it "stabilizes" the variance at small values of x provided that the transformation parameter c > 0 is properly chosen [32, 33].

In Eq. (33) the glog-transformed data are averaged over the probability distribution $p_p^P(x,r-x)$ representing the convolution product of the mutually independent probabilities of the non-specific Gaussian background and the specific expression signal (see previous section and [30]). This approach assumes that the probability distribution of each individual transcript is the same as the mean distribution averaged over all transcripts. Alternatively one can use the un-weighted kernel with p^S = const. In this case Eq. (33) restricts the averaging to the background distribution.

In the next step, the obtained values are corrected for the sequence effects according to $Z_{0,p}^P=Z_p^P-\delta A_p^{P,S}$ using the positional-dependent sensitivity model and the Z_p^P values instead of the log-intensities for parameterization (Eq. (50)). Then, the probe-level data are summarized for each probe set by calculating the Tukey-biweight median, $Z_{0,set}^P=TB{(Z_{0,p}^P)}_{p\in set}$[25], and finally by transforming them into linear scale to get the probe set-level expression degree, $L_{set}^{P,S}={10}^{Z_{set}^P}$. Eq. (33) thus provides "PMonly" and "MMonly" estimates of the expression degree with P = PM and MM, respectively.

with $\Delta L_{0,p}=\left(L_{0,p}^{PM}-L_{0,p}^{MM}\right)$ and $\Delta L_p^N(x,{\mu }_c^P)=\left(L_p^{PM,N}(x,{\mu }_c^{PM})-L_p^{MM,N}(x,{\mu }_c^{MM})\cdot {10}^{-x\cdot (1-\rho )}\right)$

Eq. (34) uses the background distribution of the PM and the conditional expectation of the bivariate normal distribution for the argument of the MM-background which gives rise to an effectively shrunken integration variable, ${\mu }_c^{MM}+x\to {\mu }_c^{MM}+{\rho }_{c}x$ (here ρ_c is the coefficient of correlation, see appendix B). Note that this approach is to some degree similar to the GC-RMA algorithm [34] which however uses pseudo-MM representing subsets of PM-probes of the same GC-content as the considered PM. A second difference to Eq. (34) is that GC-RMA directly subtracts the (non-linear) MM-intensity without explicit consideration of the non-specific background.

In this context we stress the fact that the non-specific background intensity is rather a variable contribution which progressively decreases with increasing amount of specific binding than a constant (see appendix A and the figure shown there). The background correction of the intensity (as realized, e.g. by GC-RMA) neglects this trend and therefore it is expected to over-correct the signal in the S- and as-hybridization ranges. As a consequence, this approach underestimates the expression degree due to two reasons: Firstly, this over-correction of the background and, secondly, the neglect of saturation. On the other hand, the linearized signal used here (Eq. (32)) corrects the intensity for saturation and, in addition, it contains a constant background level corrected by simple subtraction in Eqs. (33) and (34).

Note that the PM-estimate exceeds the MM signal roughly by the factor s_c ≈ 5 – 10 at comparable variance of the residual background distribution (see Figure 9, part b). For the coefficient of variation, CV ≈ σ/S, one gets roughly CV(PM) ≈ CV(PM-MM) << CV(MM). The MM consequently are expected to provide considerably less accurate expression estimates compared with the respective PMonly and PM-MM values.

The characteristic expression index (or exponent of specific binding) φ_c complements the respective exponents for non-specific hybridization β_c and the characteristic S/N-index λ_r as the basic chip summary-characteristics of specific and non-specific binding.

Chip characteristics and expression estimates in natural units

It specifies the mean occupancy of the probes in the absence of specific transcripts, Θ^PM(0) ≈ X_c^PM,N.

Hence, the hook method measures both, the abundance of specific and non-specific transcripts in the hybridization solution in chip-related units such as the relative occupancy of the probe spots and the respective binding strengths.

The specific and non-specific signals are combined into the S/N-ratio, R (Eq. (26)), which provides the expression degree in units of the non-specific background contribution in the particular chip experiment. The S/N-ratio is directly related to the hook-coordinates of a selected probe set and thus it can be roughly deduced by visual inspection of the particular hook curve (Eqs. (46) and (18)). Figure 8 illustrates the relation between the intrinsic expression measures and the hook coordinates for a typical microarray hybridization. The point of half coverage (Θ^PM = 0.5, X^{PM, S} ≈ 1) is found beyond the maximum in the sat-regime. Virtually no probe of the chosen chip meets this condition. Note that Θ^PM scales with the distance relative to the end point referring to R → ∞ (see above). Vice versa, the mean fraction of specifically hybridized oligomers, x^{PM, S}, scales with the distance relative to the N-point. It steeply increases in the mix-regime and reaches the conditions at which 50% of the bound transcripts belong to the specific ones, x^{PM, S} = 0.5, at relatively small abscissa-values. Hence, specific hybridization starts to dominate over non-specific one always at the beginning of the mix-range. The S-range near the maximum of the hook-curve refers to probes with a 50 – 100 fold excess of specific hybridization, R^PM ≈ 50 – 100. The width of the hook of about β = 2.7 is equivalent with the background strength of N_c ≈ 10^-2.7 = 0.002 which in turn rescales the S/N-ratio into binding strengths, for example S_c = 0.1–0.2 for R = 50–100. These rough estimation shows that the maximum of the hook is equivalent with the occupancy of the PM-spots in the order of 10 – 25%.

The mean binding constants, K_c^PM,S and K_c^PM,N and thus also the used measures of specific and non-specific hybridization depend on the particular probe and chip design (e.g., the length of the probe-oligonucleotides, their density and the type of the mismatches used for the MM probes). Consequently they are specific for the used chip type, on one hand. On the other hand, also the conditions of a particular hybridization affect the K_c^PM,h because their values depend on all processes affecting the binding reaction between the probe-oligonucleotides and the targets. For example, the composition of a particular sample will affect K_c^PM,N and K_c^PM,S as well, because both constants depend on the extent of target-dimerization which is a function of the concentrations of the reacting species in the hybridization solution [1].

A. Concentration scales of the Δ-vs-Σ-trajectory: Derivation of Eqs. (16) – (18)

Eq. (43) shows that the specific coverage monotonously increases with the S/N-ratio whereas the non-specific "background" coverage remains virtually constant at R < 1/X^{PM, N} for the PM (and at R < (n/(s·X^{PM, N}) for the MM, see Figure 12). Specific transcripts of large binding strength progressively replace the bound non-specific fragments with further increasing values of R. As a consequence the non-specific coverage Θ^{P, N} drops and finally disappears for R → ∞.

and Eq. (17) for the mean fraction, <x^S>. Here we make use of ${\Theta }^{\text{PM},\text{N}}(0)={10}^{-\beta -\frac{1}{2}\log \text{n}}$ and ${\Theta }^{\text{MM},\text{N}}(0)={10}^{-\beta +\frac{1}{2}\log \text{n}}$ (see Eqs. (42) and (14)).

and Eq. (18) for the respective average $〈\text{R}〉\equiv {10}^{0.5\left(\log {\text{R}}^{\text{PM}}+\log {\text{R}}^{\text{MM}}\right)}$.

B. The slope of the hook-curve in the mix- and the N-ranges

The initial slope of the mix range essentially depends on the "height"-parameter α_c and thus on the PM/MM-ratio of the specific binding constants (see Eqs. (13) and (4)).

Here σ_PM², σ_MM² and σ_PM,MM denote the variance and covariance of the set-averaged log-intensities of the PM and MM probes in the N-range. The formula for the correlation coefficient ρ assumes $<{\sigma }^2>\approx {\sigma }_{\text{PM}}^2\approx {\sigma }_{\text{MM}}^2$ to a good approximation (see Figure 9).

Hence, a tiny slope near zero indicates strongly correlated PM and MM signals (ρ~1) whereas the increase of slope(N) reflects "decoupling" of PM and MM intensities.

C. Fit of the position-dependent sensitivity profiles

Here f_k(b_m) is the probability to find the subsequence b_m at sequence position k within the probes of a probe set. The weighted least squares algorithm fits Y(theory) to Y(exp) by optimizing the coefficients δε_k^P,h(b_m) using single value decomposition [35].

D. Fit of the Langmuir model to the Δ-vs-Σ data

Equation (53) thus returns the S/N-ratio R for the Σ₀^hook-coordinate of a probe set and a given set of model-parameters {α_c, β_c, Σ_c(0), Δ_c(0)}. The parametric equation for the hook-ordinate (Eq. (10)) then provides an estimate of Δ₀(R) referring to the respective probe set. The fit-algorithm minimizes the sum of least squares between the calculated and measured Δ-values, $\text{SSQ}={\displaystyle \sum _{\text{i}}{\text{w}}_{\text{i}}\cdot {\left({\Delta }_{\text{i},0}({\text{R}}_{\text{i}})-{\Delta }_{\text{i},0}^{\text{hook}}\right)}^2}$, by adjusting the parameters α_c, β_c and Σ_c(0) using an iterative gradient method with linearly decreasing step size. The ordinate value of the start-point was set to the break between the N- and mix-range, ${\Delta }_{\text{c}}(0)={\Delta }_0^{\text{break}}$ (see above).

We divide the Σ₀^hook-axis into i = 10 – 30 equally-spaced sampling points Σ₀^hook = Σ_i by averaging over the probe/probe set data within the respective sampling interval of width δ Σ, i.e., ${\Delta }_{\text{i},0}^{\text{hook}}=〈{\Delta }_{\text{p},0}^{\text{hook}}〉$ for all probes with ${\Sigma }_{\text{i}}-\delta \Sigma \le {\Sigma }_{\text{p},0}^{\text{hook}}<{\Sigma }_{\text{i}}+\delta \Sigma$ which were used in Eq. (53) to calculate the respective S/N-ratios R_i and theoretical ordinate points Δ_i,0(R_i) using Eq. (10). The weighting factor, ${\text{w}}_{\text{i}}\approx \frac{\rho ({\Sigma }_{\text{i}})}{{\sigma }_{\Delta }{({\Sigma }_{\text{i}})}^2}$, considers the data-density, ρ(Σ_i), and variance, σ_Δ(Σ_i)², per interval. Both, ρ(Σ_i) and σ_Δ(Σ_i)², are given by decaying functions ([30, 32]) which partly compensate each other to a good approximation. By default the weighting factor is therefore set to w_i = 1.