An efficient algorithm for the stochastic simulation of the hybridization of DNA to microarrays

Abstract

Background: Although oligonucleotide microarray technology is ubiquitous in genomic research, reproducibility and standardization of expression measurements still concern many researchers. Cross-hybridization between microarray probes and non-target ssDNA has been implicated as a primary factor in sensitivity and selectivity loss. Since hybridization is a chemical process, it may be modeled at a population-level using a combination of material balance equations and thermodynamics. However, the hybridization reaction network may be exceptionally large for commercial arrays, which often possess at least one reporter per transcript. Quantification of the kinetics and equilibrium of exceptionally large chemical systems of this type is numerically infeasible with customary approaches.

Results: In this paper, we present a robust and computationally efficient algorithm for the simulation of hybridization processes underlying microarray assays. Our method may be utilized to identify the extent to which nucleic acid targets (e.g. cDNA) will cross-hybridize with probes, and by extension, characterize probe robustnessusing the information specified by MAGE-TAB. Using this algorithm, we characterize cross-hybridization in a modified commercial microarray assay.

Conclusions: By integrating stochastic simulation with thermodynamic prediction tools for DNA hybridization, one may robustly and rapidly characterize of the selectivity of a proposed microarray design at the probe and "system" levels. Our code is available at http://www.laurenzi.net.

Figure 1

Structure of the hybridization table (HT). Rate constants and populations are stored in such a way as to simplify the calculation of the quantities α_j, ϕ_j, α and Φ (See Eqs. 7, 8, 9, and 10) thereby reducing the number of operations per time step of the simulation from O(2N_P N_T) to O(N_P + N_T).

Figure 2

Determination of equilibrium within simulations of hybridization. Results are shown for the hybridization of all 6256 Agilent probes with 6718 full-length yeast cDNAs. The red circle indicates the time at which equilibrium is attained.

Figure 3

Effect of kinetic rate constants upon the equilibrium state in microarray simulations. The fractional occupancies of probes at equilibrium (Eq. 27) are unaffected by random perturbations of rate constants as t → ∞. Each point represents a pair of results for each of the 6256 Agilent probes targeting yeast ORFs. The CPU time required for equilibration of simulated systems (4 h, 100 h, above), like the hybridization time (not shown) does, however, depend upon the rates. Our methodology for estimating rate constants yields rapidly converging simulations.

Figure 4

Comparison of SSAs and the Law of Mass Action. Results of simulations of the hybridization of ten Agilent probes with their targets (N_P = N_T = 10) under standard experimental conditions are illustrated (mean ± SD, n = 5). All three approaches yield statistically indistinguishable results. In these simulations there are 1000 probe molecules per feature (N = 1000), corresponding to a hybridization volume of 0.275 nL.

Figure 5

Computational performance of the Method of Partial Sums and Next Reaction Method. CPU times for of simulations of hybridization until equilibrium are illustrated. The in silico time t required to establish equilibrium is determined by our method and utilized as an end point in "Next Reaction" simulations. Our algorithm outperforms the Next Reaction Method in both absolute terms as well as on a per-probe basis (frame).

Figure 6

Hybridization of cDNA with Agilent 4 × 44 array probes. The populations of hybrids with these fourteen probes (mean ± SD) are a subset of the complete set of results, which are generated from five replicate simulations.

Figure 7

Selectivity distribution for the Agilent probe set. The black line illustrated the selectivities (Eq. eq:Selectivity) of all Agilent probes when the array is hybridized to yeast cDNAs at the initial concentrations described in Methods. The gray lines are results for four different initial initial conditions. The red line delimits 90% selectivity. The distribution is independent of the initial concentrations (x-axes are different for all five sets of results). About 100 probes will not be selective in any given microarray experiment.