BGX: a Bioconductor package for the Bayesian integrated analysis of Affymetrix GeneChips
 Ernest Turro^{1}Email author,
 Natalia Bochkina^{2},
 AnneMette K Hein^{3} and
 Sylvia Richardson^{1}
DOI: 10.1186/147121058439
© Turro et al; licensee BioMed Central Ltd. 2007
Received: 11 September 2007
Accepted: 12 November 2007
Published: 12 November 2007
Abstract
Background
Affymetrix 3' GeneChip microarrays are widely used to profile the expression of thousands of genes simultaneously. They differ from many other microarray types in that GeneChips are hybridised using a single labelled extract and because they contain multiple 'match' and 'mismatch' sequences for each transcript. Most algorithms extract the signal from GeneChip experiments in a sequence of separate steps, including background correction and normalisation, which inhibits the simultaneous use of all available information. They principally provide a point estimate of gene expression and, in contrast to BGX, do not fully integrate the uncertainty arising from potentially heterogeneous responses of the probes.
Results
BGX is a new Bioconductor R package that implements an integrated Bayesian approach to the analysis of 3' GeneChip data. The software takes into account additive and multiplicative error, nonspecific hybridisation and replicate summarisation in the spirit of the model outlined in [1]. It also provides a posterior distribution for the expression of each gene. Moreover, BGX can take into account probe affinity effects from probe sequence information where available. The package employs a novel adaptive Markov chain Monte Carlo (MCMC) algorithm that raises considerably the efficiency with which the posterior distributions are sampled from. Finally, BGX incorporates various ways to analyse the results, such as ranking genes by expression level as well as statistically based methods for estimating the amount of up and down regulated genes between two conditions.
Conclusion
BGX performs well relative to other widely used methods at estimating expression levels and fold changes. It has the advantage that it provides a statistically sound measure of uncertainty for its estimates. BGX includes various analysis functions to visualise and exploit the rich output that is produced by the Bayesian model.
Background
Oligonucleotide microarrays allow biomedical researchers to estimate the expression of thousands of genes simultaneously through their mRNA transcripts. A labelled, fragmented version of the RNA may be hybridised onto an array containing hundreds of thousands of complementary oligonucleotides and then scanned. Affymetrix 3' GeneChip arrays represent genes by sets of probe pairs, each of which consists of an oligonucleotide of length 25 which matches a corresponding RNA subsequence perfectly (PM) and an identical probe with an inverted oligonucleotide on position 13 (MM) that is intended to measure nonspecific hybridisation.
The BGX model [1] is an integrated approach to the analysis of GeneChip microarrays in which correction for nonspecific hybridisation and gene expression level estimation are performed simultaneously. Posterior distributions of parameters in the model may be obtained numerically. Based on these distributions, a powerful method for detecting differential expression has been developed [2].
The probes on Affymetrix GeneChips have been found to exhibit varying propensities to "shine" according to the base composition of their sequences [3] and methods for estimating expression levels from GeneChips that incorporate probe affinity effects have shown demonstrable advances over methods in which these effects are ignored (see, e.g. [4]). We present a new Bioconductor [5] package that implements the BGX model, includes an extension to incorporate probe affinity effects, employs novel algorithmic techniques to sample effectively from posterior distributions, and provides various analysis and plotting functions.
Implementation
Basic model
where TN denotes the truncated normal distribution, truncated to the positive axis. The central parameter of equation (3), μ_{ gc }, acts as the BGX expression measure, and equations (1) to (4) represent the basic BGX model.
The core of the model is implemented in the C++ programming language for efficiency and uses MCMC to sample from the full posterior distributions of each parameter. Parameters are estimated using Gibbs sampling where possible (φ and τ) and a Random Walk MetropolisHastings algorithm elsewhere (S, H, μ, σ, λ and η). Three C++ class templates are used to instantiate zero, one and twodimensional MCMC update objects for each parameter according to the dimensionality of the corresponding suffixes. Each of the instantiated objects is updated in sequence using references to all other necessary parameters during a burnin period, which is discarded, and a sampling period, which is used for the posterior distributions.
Probe affinity extension
The extended model, which we denote GCBGX, is based on equations (1), (2), (3) and (5).
The probes on the arrays are, prior to analysis, grouped into a number of probe affinity categories. This is done by: (a) calculating the probe affinities using the gcrma Bioconductor package [4], (b) rounding them to the first decimal place, (c) assigning each value to a preliminary probe affinity category and (d) ensuring that the final categories contain a sufficient number of probes by collapsing small preliminary categories together. We enumerate the resulting probe affinity categories 1, ... , K by increasing affinity. Once categorised, the affinityspecific parameters are estimated from the data, simultaneously with all other parameters.
In some cases, Affymetrix do not directly provide the sequences for all probesets due to licensing restrictions and, consequently, there are Bioconductor probe packages that do not contain complete sequence information. For example, hgu95aprobe (version ≤ 1.16.2) lacks sequences for probes belonging to 172 probesets. We tackled this problem by treating α(g, j) as a random variable, taking values from 1 to K, with prior probability equal to the observed frequency of the categories, ${p}_{k}=\frac{{N}_{k}}{N}$, where N_{ k }is the number of probes in category k and $N={\displaystyle {\sum}_{k=1}^{K}{N}_{k}}$.
Adaptive MCMC
The full conditional distributions of S, H, μ, σ, λ and η are updated by drawing new values from a proposal distribution, typically a Random Walk (RW) Gaussian proposal centred on the current value with a chosen variance. A typical experiment consists of several hundred thousand probes, resulting in potentially millions of S and H components and tens of thousands of μ and σ components. Each component of a given parameter has a different support and consequently a different optimal RW proposal variance. Using a fixed variance for all components results in excessively low or high acceptance ratios for a large proportion of components, leading to highly autocorrelated chains.
In order to tackle this problem, we implemented the novel Adaptive MetropolisWithinGibbs algorithm recently proposed by Roberts and Rosenthal [6, 7]. We used a unique proposal variance for each object, which adapts to its optimal value after successive batches of 50 iterations. The aim is to achieve an acceptance ratio of around 0.44, which has been shown to be optimal for onedimensional proposals in certain settings [8, 9], and is commonly accepted as being a sensible benchmark. An acceptance rate that is close to zero implies inefficient mixing, while an acceptance rate that is close to one implies the probability space is not efficiently explored. The algorithm proceeds as follows:

For each component c of parameter p, assign a parameterspecific starting value to the corresponding proposal variance, ${\sigma}_{cp}^{2}$

Choose a sequence δ (n) → 0. We chose δ (n) = min(0.01, n^{1/2})

Start the MCMC simulation

After the n^{th} batch of 50 iterations, calculate the acceptance ratio over the last batch

If the acceptance ratio is less than the optimal value of 0.44, increase log(${\sigma}_{cp}^{2}$) by δ (n), else decrease it by δ (n)
The algorithm preserves ergodicity as long as each kernel has the right stationary distribution; the total variation distance between successive kernels tends to zero in probability; and the convergence time of each kernel is bounded in probability [6].
R package
The C++ component of BGX is compiled as a shared object which is loaded and executed automatically from within the R package [10]. BGX integrates standard Bioconductor classes such as AffyBatch to store raw microarray data and ExpressionSet to store processed gene expression measures. Users interested in running BGX programmatically from a shell script, for instance, or in a more memoryefficient manner, also have the choice to run a standalone binary version of the program.
Results and Discussion
Usage
The bgx package and its dependencies, affy and gcrma, may be installed automatically from the Bioconductor repository from an R shell. The package contains documentation and executable examples in a "vignette" file available using openVignette(). Users who wish to compile bgx from source will require the Boost C++ libraries [11] and the hgu95av2cdf Bioconductor package. The core functionality of the package is contained in the bgx function, which takes an AffyBatch object instantiated from one or more GeneChip CEL files as its first argument and returns an ExpressionSet object containing expression values for each gene and condition:
aData < ReadAffy("chip1.CEL","chip2.CEL")
eset < bgx(aData)
assayData(eset)$exprs # Returns expression values
assayData(eset)$se.exprs # Returns standard errors for expression values
Optional arguments include samplesets, which specifies the experimental design; genes, which specifies a subset of genes to analyse; burnin and iter, which specify the number of iterations for the burnin and post burnin phases of the algorithm respectively; probeAff, which specifies whether or not to use the probe affinity extension to the original BGX model, and adaptive, which specifies whether or not to use MetropolisWithinGibbs step adaptation. Full documentation for the bgx function is available by running help(bgx).
Although the point measures returned in the ExpressionSet object are useful, the distinctive power of the BGX method is that it provides samples from the full posterior distributions of the expression parameter, μ_{ gc }. These samples are, by default, saved in directories named run.1, run.2, etc. in R's current working directory, although this may be overridden with the rundir argument. They may be read into R in order to analyse the results of a simulation as follows:
bgxOutput < readOutput.bgx("run.1")
For the purposes of this paper, BGX was run with the "gold standard" of 16 k burnin iterations and 64 k sampling iterations. However, the recommended 8 k burnin iterations and 16 k sampling iterations are sufficient to provide good estimates of μ_{ g }(Additional File 1). Under these settings, BGX takes approximately one hour per array on a standard 64bit 3 GHz computer. Analyses of up to 100 arrays ought to "fit" in a computer equipped with 4 GB of memory. However, BGX may be run separately for each condition, and the output subsequently combined in R by passing multiple output directories to the readOutput.bgx function. Since φ, the only parameter that is shared between conditions, is very stable for a given type of array, the impact on the output of running BGX separately on each condition is negligible.
Estimation of nonspecific hybridisation
Performance of adaptive MCMC
Evidently, if the sample is not autocorrelated, then P_{ j }= 0 for all j, a(X) collapses to 1 and var($\overline{X}$) becomes equal to the familiar expression for an IID sample, E [S^{2}(X)/n]. From (6), the IACT of a chain relates positively with the variability of its mean, and thus highly autocorrelated chains lead to poor estimates of our gene expression measure.
One way of improving our estimates is to increase the number of iterations while maintaining a fixed subsample size. This translates to subsamples being further apart on the original chain and therefore less correlated. It is faster and more attractive, however, to use an adaptive algorithm that explores the probability space more efficiently. Using the Golden Spike data set [12] for our investigation, we found that the adaptive method led to a range of optimal proposal magnitudes for the MetropolisHastings parameters. Figure 5 illustrates this with a histogram of the optimal log variance for S proposals on one array and the original fixed step size overlaid in black. Figure 6 (left & centre) shows a dramatic reduction in the IACT of the S parameters and a milder improvement on the μ_{ g }parameters of expressed genes. A similar improvement was observed for the IACT of the H parameters, this time for all genes (Figure 6 right).
dg = μ_{g 2} μ_{g 1}, g = 1 ..., G, that is, samples from the posterior distribution of the difference in the BGX expression measure for each gene, and $\stackrel{\_}{a({d}_{g})}$, the estimate of the Monte Carlo standard error, is calculated using Sokal's adaptive truncated periodogram estimator [?]. Our zscore differs from the measure used in [2] by a factor of $\sqrt{((n/\stackrel{\_}{a({d}_{g})}1)/(n1)}$, which takes into account the autocorrelation structure of the sequence of values generated by the algorithm. Since the adaptive MCMC algorithm has the effect of decreasing a(d_{ g }) for expressed genes while keeping it approximately constant for nonexpressed genes, it leads to an increase in the ranking of expressed genes and consequently in BGX's capacity to detect differential expression.
Performance on spikein datasets
We illustrate the performance of bgx by presenting detailed results from analyses of arrays from the Affymetrix Latin Square data [15] and the Golden Spike data set [12].
Latin Square data
Affymetrix published two data sets for assessing the performance of expression algorithms on their microarrays. The HGU95A data set consists of 16 genes spiked in at known concentrations ranging from 0 to 1024 pM and arrayed in a Latin Square format. We considered 16 instead of the original 14 genes described by Affymetrix because we included two extra spikeins, 546_at and 33818_at, as reported in [16]. We used two replicates and 14 unique concentration configurations labelled A to M and Q. 2716 of the probes in this data set had no sequence information and therefore their probe affinity categories were estimated from the data as part of the model. The HGU133A data set consists of 64 genes spiked in at known concentrations ranging from 0 to 512 pM. We considered 64 instead of the original 42 genes described by Affymetrix because we included 22 extra spikeins, as reported in [17]. We used all 3 replicates for each of the 14 concentration groups.
Golden Spike data set
The Golden Spike data set consists of six DrosGenome1 GeneChips, with three technical replicates from two conditions: C and S. There are 14010 probe sets in each array representing 14010 genes. 2535 of these are expressed equally under both conditions while 1331 genes are upregulated in S relative to C. The data is highly valuable for comparing chip analysis methods because it is fully controlled and contains very realistic noise. Due to the asymmetry of the spikeins, a normalisation of the posterior distributions similar to that advocated in [12] was carried out by fitting a loess curve to the MA plot [19] of the posterior mean values of μ_{ gc }for the nondifferentially expressed genes, predicting a curve from the fit for all genes, and subtracting the curve from the posterior distributions of the differences in expression. The RMA, GCRMA and MAS5 expression measures were similarly adjusted using loess normalisation at the probeset level instead of the default quantile normalisation at the probe level.
Conclusion
BGX is a new Bioconductor R package for analysing 3' Affymetrix GeneChips. BGX implements a fully integrated Bayesian hierarchical model with the option to take into account sequencedependent probe affinities. BGX uses a novel adaptive MCMC algorithm that improves the efficiency with which the posterior distributions of parameters are sampled from. BGX compares favourably to RMA and GCRMA at detecting differential expression, particularly at low concentration levels.
Availability and requirements
Project name: BGX
Project homepage: http://bgx.org.uk
Operating systems: Platform independent
Programming language: C++, R
Other requirements: R, Bioconductor
License: GNU GPL
Any restrictions to use by nonacademics: No
Declarations
Acknowledgements
We thank Jeffrey S. Rosenthal, Gareth O. Roberts for advice on adaptive MetropolisHastings algorithms, and Alex Lewin and Marta Blangiardo for valuable discussion. This work was supported by BBSRC 'Exploiting Genomics' grant 28EGM16093 and the John and Birthe Meyer Foundation.
Authors’ Affiliations
References
 Hein AMK, Richardson S, Causton HC, Ambler GK, Green PJ: BGX: a fully Bayesian integrated approach to the analysis of Affymetrix GeneChip data. Biostatistics 2005, 6(3):349–373. 10.1093/biostatistics/kxi016View ArticlePubMedGoogle Scholar
 Hein AMK, Richardson S: A powerful method for detecting differentially expressed genes from GeneChip arrays that does not require replicates. BMC Bioinformatics 2006, 7: 353. 10.1186/147121057353PubMed CentralView ArticlePubMedGoogle Scholar
 Naef F, Magnasco MO: Solving the riddle of the bright mismatches: labeling and effective binding in oligonucleotide arrays. Phys Rev E Stat Nonlin Soft Matter Phys 2003, 68(1 Pt 1):011906.View ArticlePubMedGoogle Scholar
 Wu Z, Irizarry RA, Gentleman R, MartinezMurillo F, Spencer F: A ModelBased Background Adjustment for Oligonucleotide Expression Arrays. Journal of the American Statistical Association 2004, 99: 909–917(9). 10.1198/016214504000000683View ArticleGoogle Scholar
 Bioconductor[http://bioconductor.org]
 Roberts GO, Rosenthal JS: Coupling and Ergodicity of Adaptive MCMC. Journal of Applied Probability 2005. [To appear] [http://www.probability.ca/jeff/ftpdir/adapt.pdf] [To appear]Google Scholar
 Roberts GO, Rosenthal JS: Examples of Adaptive MCMC.2006. [http://probability.ca/jeff/ftpdir/adaptex.pdf]Google Scholar
 Roberts GO, Gelman A, Gilks WR: Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms. The Annals of Applied Probability 1997, 7: 110–120. 10.1214/aoap/1034625254View ArticleGoogle Scholar
 Roberts GO, Rosenthal JS: Optimal Scaling for Various MetropolisHastings Algorithms. Statistical Science 2001, 16(4):351–367. 10.1214/ss/1015346320View ArticleGoogle Scholar
 The R Project for Statistical Computing[http://www.rproject.org]
 Boost C++ Libraries[http://boost.org]
 Choe SE, Boutros M, Michelson AM, Church GM, Halfon MS: Preferred analysis methods for Affymetrix GeneChips revealed by a wholly defined control dataset. Genome Biol 2005, 6(2):R16. 10.1186/gb200562r16PubMed CentralView ArticlePubMedGoogle Scholar
 Anderson T: The Statistical Analysis of Time Series. New York: John Wiley and Sons; 1971.Google Scholar
 Law AM: Confidence intervals in discrete event simulation: a comparison of replication and batch means. Naval Res Logistics Quart 1977, 23: 667–678. 10.1002/nav.3800240414View ArticleGoogle Scholar
 Affymetrix – Latin Square Data[http://www.affymetrix.com/support/technical/sample_data/datasets.affx]
 Cope LM, Irizarry RA, Jaffee HA, Wu Z, Speed TP: A benchmark for Affymetrix GeneChip expression measures. Bioinformatics 2004, 20(3):323–331. 10.1093/bioinformatics/btg410View ArticlePubMedGoogle Scholar
 McGee M, Chen Z: New SpikedIn Probe Sets for the Affymetrix HGU133A Latin Square Experiment. COBRA Preprint Series 2006.Google Scholar
 Irizarry RA, Hobbs B, Collin F, BeazerBarclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 2003, 4(2):249–264. 10.1093/biostatistics/4.2.249View ArticlePubMedGoogle Scholar
 Yang YH, Dudoit S, Luu P, Lin DM, Peng V, Ngai J, Speed TP: Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nucl Acids Res 2002, 30(4):e15. 10.1093/nar/30.4.e15PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.