BACKGROUND
affymetrix genechips are one of the most widely used commercially available oligonucleotide arrays.
they have gained widespread popularity for a number of reasons, among which are their high degree of standardization of the production process and their ability to interrogate tens of thousands of genes simultaneously.
they differ from many other array types in that they are one color  arrays and because each gene is represented by a probe set of pairs of perfect match  and mismatch  probes.
each pm is chosen to match a particular  <dig> base pair stretch of the sequence encoding the gene, and the full set of pms is chosen with the aim of uniquely identifying the gene.
the accompanying mms are identical to their pm counterparts except for a complementary base substitution at the middle nucleotide.
they are intended to be used to correct for non-specific hybridization.
the full set of pms and mms for a gene represent the basis for the estimation of the level of expression of the gene.
most genechip based expression studies are carried out with few replicates.
there are two major factors behind this: the considerable cost of the genechip arrays and limitations on the amount of available rna.
as microarray expression studies are prone to experimental imperfections, failed hybridizations are often encountered resulting in further reduction in the number of replicates available for analysis.
if the number of replicates falls below three, most available analysis tools become unsuitable or unapplicable because they rely on the estimation of variances which is difficult in such circumstances.
thus, the development of methods for analyzing experiments with few or no replicates is of high importance.
bayesian gene expression, bgx  <cit> , is an integrated approach to the analysis of affymetrix genechip arrays.
it relies on the formulation of a bayesian hierarchical model for estimating expression levels from probe level genechip data.
in the bgx approach background correction, gene expression level estimation and differential expression is performed in an integrated analysis, allowing all uncertainties to be dealt with simultaneously in a coherent statistical framework.
posterior distributions of gene and condition specific bgx expression levels are obtained from a simultaneous consideration of the probe pair intensities in the available probe sets representing the gene.
samples from the posterior distributions are generated using markov chain monte carlo methods.
if replicate arrays are available the information in their probe set intensities will be considered jointly in the estimation of expression levels.
replicate arrays, however, are not essential for obtaining the posterior distributions of expression for the genes – these will be obtained from the collection of intensities on the array even when only a single array is available.
samples from the posterior distribution of the differences in expression levels present a basis for inference on differential expression.
in this paper we develop a method for performing differential expression studies from genechip experiments without replicates.
the procedure exploits the posterior distributions of differences in expression, obtained from a bgx analysis of genechip arrays, for creating ranked gene lists.
in order to define suitable cut-offs for the list, an estimate of the proportion of differentially expressed genes is obtained by empirically estimating the null distribution of a relevant statistic using an approach similar to that of efron  <cit> .
the performance of the method is tested on publicly available spike-in data sets and compared to those of other methods, and further evaluated in a biological study.
RESULTS
the bgx model and methodology
most methods for analyzing genechip arrays adopt a stepwise procedure for obtaining a point estimate of expression for each gene on each array.
the steps in the procedures typically consist of background correction and normalization followed by summarization  or the fitting of a linear model, often performed on the log-scale background corrected intensities .
having obtained a point estimate of expression for each gene on each array, studies of differential expression between pairs of conditions are carried out by comparing the collections of point estimates under the conditions using a t-type statistic such as in sam  <cit> , cyber-t  <cit>  or limma  <cit> .
the bgx method differs from these stepwise point estimate approaches in that  all steps in the analysis are dealt with simultaneously,  gene and condition specific expression levels are estimated from a joint consideration of the available probe set intensities and  the outcomes are posterior distributions of expression rather than point estimates.
thus, uncertainties associated with each of the steps are taken into account at all levels of analysis, and the joint uncertainty on the expression level for a gene under a condition is reflected in the shape of the posterior distribution obtained for the level of expression for that gene under that condition.
explicitly, with g, j, c and r denoting gene, probe, condition and replicate, respectively, let sgjcr and hgjcr denote gene, probe, condition and replicate specific and non-specific binding  of rna, and let φ ∈  be a fraction.
to further allow for additive array specific errors, e.g.
accommodating mms bigger than pms, the bgx model hypothesizes:
pmgjcr ~ n
mmgjcr ~ n.
information on the level of expression of gene g under condition c is represented by the set of signal parameters representing the gene under this condition: sgjcr, j =  <dig> ...,j, r =  <dig> ..., rc.
we assume that these, shifted and logged, come from a gene and condition specific truncated normal distribution.
the non-specific hybridization parameters hgjcr reflect characteristics specific to the sample hybridized, leading us to assume array specific truncated normal distributions for these, shifted and logged.
thus,
log ~ tn,
log ~ tn.
we will here refer to the μgc parameters as the bgx expression indices or levels.
we assume exchangeability of the gene and condition specific variance parameters,
log() ~ n,
with ac and  fixed at values obtained by an empirical bayes like approach, thus stabilizing the variance estimation.
in all of the above, the distributions are conditional on variables on the right hand side, and independent for all suffices.
the model is fully specified by declaring the following, generally weakly informative priors, independent for all suffices: μgc ~ u, φ ~ , λgc ~ n),  ~ Γ and  <dig> ~ Γ).
for a more in-depth discussion of the model we refer to  <cit> .
the bgx model relies on mcmc methods for obtaining samples from the posterior distributions of the parameters.
the shapes of the posterior distributions of the bgx gene expression indices, μgc, are determined by the probe response patterns .
thus, a highly consistent probe set response leads to a tight posterior distribution of expression, and a less consistent pattern will result in a flatter, possibly multi-modal, posterior distribution.
examples of posterior expression index distributions, μg, <dig> and μg, <dig>  are given in figure  <dig> .
the corresponding kernel density plots for the differences in expression indices, dg = μg, <dig> - μg, <dig>  are given in figure  <dig> .
the uncertainties of the expression indices are reflected in the shape of these distributions.
for gene  <dig> the multi-modality of the posterior distribution of the expression index under condition  <dig>  is reiterated by the multi-modal posterior distribution of the difference in expression.
for the other two example genes the posterior distributions of the differences in expression are tight and uni-modal, centered close to zero and around one respectively, indicating similar expression levels for gene  <dig> and different expression levels for gene  <dig> under the two conditions considered.
addressing differential expression with replicates
a popular approach to conducting differential gene expression studies is to rank the genes according to their degree of evidence for differential expression, and to estimate false discovery rates for different cut-offs on the ranked gene list.
this allows the experimenter to obtain a prioritized list of genes to pursue in follow-up studies, with a guidance as to how many genes on the list are expected to be false positives.
such approaches are taken in the implementations of the sam, limma and cyber-t methods.
each of the methods calculate a different modified t-statistic, the modification relating to the standard deviation or variance calculation in the denominator, and genes are ranked on the resulting p-values.
in sam a false discovery rate is estimated based on permuting the original data to get the distribution of  t-statistics under the null-hypothesis of no differential expression.
limma is implemented with the benjamini and hochberg method  <cit>  for estimating fdr and calculation of adjusted p-values.
cyber-t adopts the method of allison et al.
<cit>  for fitting a mixture of beta-distributions  distribution) to the observed p-values, and reports estimated true and false positives along with the posterior probability of differential expression.
thus, all methods make use of point estimates of expression and depend upon replicates being available for estimation of the variance in the modified t-statistics .
addressing differential expression without replicates
without replicates the above methods for analysis of differential expression are unapplicable and alternative methods are needed.
in this section we describe how the bgx model and methodology may be exploited to obtain such a procedure.
we use features of the samples from the posterior distributions collected in the mcmc sampling to produce ranked gene lists.
the ranking takes into account the estimated difference in expression level as well as the associated uncertainty.
we then consider the set of posterior probabilities of expression differences being smaller than zero, p, g =  <dig> ..., g. by comparing their observed distribution to that expected under the null-hypothesis of no differential expression, we obtain an estimate of the number of differentially expressed genes, gd.
this allows us to choose the cut-off of the ranked gene list in an informed manner.
the procedures for ranking the genes and for estimating gd in the bgx framework are described in two separate subsections below.
the final part of this section contains a comparison of the performance of the bgx based method for performing differential expression analysis from genechips without replicates to those of other available methods: the ebarrays method of kendziorski et al.
<cit> , the wilcoxon signed rank test for comparison calls of the microarraysuite software  <cit> , and the efron  <cit>  method using the standardized bgx differences  as z-statistics.
ranking genes using bgx
in the bgx framework, the samples from the posterior distributions of the differences dg = μg, <dig> - μg, <dig>  g =  <dig> ...,g, represent a natural base for inference on differential expression between conditions  <dig> and  <dig>  these are available irrespectively of the number of replicates for each condition used in the analysis.
there are numerous ways in which these posterior distributions can be exploited with the aim of addressing expression differences.
here we study two types of rankings reflecting the potential of the genes as promising candidates for differential expression:  ranking on the 'standardized bgx differences', zg = mean/sd, where the mean and standard deviation are computed from the posterior sample of dg values, and  ranking on the highest percentile, α*, for which the α-percent credibility interval for the difference dg does not cover zero.
note that both rankings use the levels of differential expression  as well as the uncertainty of these  in the ranking.
without replicates point estimate based methods clearly do not have this ability.
to illustrate ranking , consider the posterior distributions of expression index differences in figure  <dig> .
all sampled values from the posterior distribution of d <dig> are above zero and α* for gene  <dig> is indistinguishable from 100%.
for gene  <dig>  only very tight credibility intervals exclude zero, and α* for this gene is small .
for gene  <dig>  α* is around 75%.
rankings of type  and  differ in their emphasis of the posterior distribution characteristics: type  summarizes the full distribution by a traditional t-type statistic  and we would expect it to perform well when ranking gaussian-shape posterior dg distributions  and possibly less well in the presence of asymmetric or multi-modal distributions .
type  rankings use the tails of the posterior dg distributions and should deal equally well with symmetric and asymmetric distributions.
however, we rely on a finite sample from the mcmc procedure  for approximating the distributions, and the estimation of the tails of the distributions may be fragile.
we examined the performance of the two ranking procedures for the bgx method on the results from the nine analyses of pairs of arrays from the choe data set  <cit> , that consist of a condition c array and a condition s array.
on each pair of arrays we ran the bgx model with  <dig> conditions and  <dig> replicate per condition, and obtained ranked gene lists using each of the rankings  and .
for the choe data it is known exactly which genes are differentially expressed .
this allows the exact numbers of true and false positives and negatives to be calculated, for all possible cut-offs in ranks of the ranked gene lists.
the two rankings resulted in almost identical counts, and there is no indication from this analysis that either is to be preferred.
roc curves summarizing the counts obtained with ranking  are shown in figure  <dig>
the performance of the bgx-based rankings from the single replicate comparisons in the choe data set is remarkable: a quarter of the  <dig> truly differentially expressed genes are included in gene lists with realized false discovery rates of  <dig>  .
by extending the gene list to  <dig> genes  the proportion of truly differentially expressed genes detected is increased to 50% and the realized false discovery rate thus about 6%.
gene lists of lengths  <dig> include 70% of the truly differentially expressed genes and have observed false discovery rates of approximately 30%.
furthermore, the curves for the nine different analyses of pairs of arrays are highly similar, indicating a stable performance.
for comparison, roc curves obtained for the same pairwise analyses using mas <dig> and rma are also given in figure  <dig>  for these analyses, genes were ranked on the absolute values of the differences in expression measures  between conditions.
with just one replicate array per condition, uncertainty of the estimates of expression is not accounted for by these methods, and they both do less well than the bgx-based method.
due to the high variability of mas <dig>  this method performs particularly poorly in the one-versus-one array rankings.
to put the above performances into perspective the roc curve obtained for the same method of ranking but using results from a bgx analysis that uses all the available arrays  are also shown in figure  <dig>  for this analysis the number of truly differentially expressed genes included in the lists are increased by  <dig>   <dig> and 5% respectively to  <dig>   <dig> and 75% for the levels of realized false discovery rates of  <dig> ,  <dig>  and  <dig> .
thus, as expected, gene ranking is improved when all replicates are used.
however, the proportion of the information contained in the full data set, that may be extracted from the single replicate analyses array is considerable.
we note that the performance of the bgx multiple array model on the full choe data set is among the methods found to perform best .
results from a similar analysis on the much less extensive affyu133a spike-in data set  <cit> , are shown in figure  <dig>  as only  <dig> of the  <dig> genes represented in this data set are differentially expressed, we plot absolute rather than fractional values of true and false positives.
for this data set the ranking of the top genes produced by the three methods are similar.
for gene lists longer than  <dig>  bgx and rma outperform mas <dig>  and rma performs better than bgx for lengths above  <dig>  thus, the relative performances of the methods differ for the choe and the affyu133a data sets.
this is most likely due to the different levels of noise in the two data sets.
with only  <dig> spike-in genes in the affyu133a data set arrays from two different conditions are almost like technical replicates.
in the choe data set, the  <dig> genes spiked-in at varying concentrations at small to moderate fold changes result in a noticeable, and more biologically realistic, level of noise between the two conditions compared.
thus, for the choe data set, noise has an impact and accounting for this, as is done with the bgx based ranking, is essential, whereas for the affyu133a data set, the impact of noise is negligible and the importance of taking this into account is outweighed by that of reproducibility of point measures for  technically replicate arrays.
estimation of the proportion of differentially expressed genes
having obtained a ranked list of genes the next question is whether we can choose a suitable cut-off.
depending on the downstream goal, we may wish to arrive at a  list that has a good chance of containing most or all of the "interesting" genes , accepting that false positives will also be included, or we may prefer to end up with a  list of genes most of which would appear to be promising candidates for differential expression, expecting very few false positives to be included.
to guide the choice of cut-off, it is useful to obtain an estimate of the proportion of differentially expressed genes.
to do this, we estimate empirically the distribution of a relevant statistic under the null hypothesis, which in turn allows for the quantification of the proportion of non-null behaving genes, following the idea of efron  <cit> .
we consider, for each gene g, the posterior distribution of the difference in expression, dg = μg, <dig> - μg, <dig>  obtained from the bgx analysis.
under the null-hypothesis of no differential expression this posterior distribution should be centered on zero.
rephrasing in terms of the posterior p probability, under the null we expect this to be  <dig> .
considering the posterior p probabilities for the full set of genes analyzed, under the null hypothesis of no differential expression there should be decreasingly less support for values away from  <dig>  towards  <dig> and  <dig>  thus, with no differentially expressed genes a histogram of the p probabilities should be uni-modal with mode of approximately  <dig> , and have smoothly decreasing tails.
the width of the central component of the histogram will depend on the posterior distributions of the dg's: the less clean the information on the expression levels , the less tight the dg distributions, and thus the flatter the histogram of p values.
values but could of course equally well have phrased them in terms of the p values).
deviations in a histogram of p values from the expected shape under the null hypothesis of no differential expression indicate the presence of differentially expressed genes: excess of p values near zero and one indicate over-expressed and under-expressed genes in condition  <dig> relative to condition  <dig>  respectively.
to quantify the number of non-null genes, we adopt an approach similar to that of efron  <cit> .
we fit a polynomial, f, to the histogram counts by poisson regression and use the central part of the histogram to estimate the null component, while excess area in the tails will represent differentially expressed genes.
to be precise, we use the following procedure : we identify the inner global max and left and right most minima on the fitted curve f, and refer to the histogram categories of these points as cmax,  and  respectively.
we use the central part of the histogram  to obtain the distribution of the p values under the null: we fit two new curves,  and , to the parts of the histogram left and right of category cmax, respectively, using the same fitting procedure as for the curve f but fixing the histogram counts of the outermost categories c <dig> and ck to zero , and giving zero weight to categories c <dig> ...,  and ,..,ck- <dig> respectively.
an estimate of the number of genes under the null is obtained by summing the fitted values of the empirical null distribution f <dig>
denoting the total number of genes by g, the estimated proportion of differentially expressed genes is
estimated numbers of over- and under-expressed genes,  and  , are obtained by quantifying the excess genes in either tail of the histogram relative to those expected under the null.
denoting the histogram count in category ci by h we set
we applied the above method to the bgx analyses of the pairs ofarrays from the choe and affyu133a data sets.
an example of a histogram of p values from a within condition comparison of two arrays from the choe data set is given in figure  <dig>  upper panel.
the full set of histograms for within condition analyses of pairs of choe arrays are summarized in figure  <dig>  left, in terms of plots of the curves, f, fitted to the histogram counts by poisson regression.
for within condition comparisons the arrays are replicates, so there should be no differentially expressed genes, and the plots indeed exhibit the shape expected under the null hypothesis: they have a central mode near  <dig>  and smoothly decreasing tails.
they thus confirm our expectations and indicate that the method works well under the null.
figure  <dig>  lower panel, and figure  <dig>  right panel, display the equivalent plots for pairwise between condition analyses of the choe data set arrays.
the histograms for these analyses exhibit a clear deviation from the shape expected under the null in terms of an excess of small p values near zero indicating the presence of over-expressed genes.
the right-hand tails of the histograms decrease smoothly and there is no indication of under-expressed genes.
this is exactly the pattern that should emerge for the choe between condition comparisons: all differentially spiked-in genes have higher concentration under condition s  than under condition c .
the estimates of the numbers of differentially expressed genes obtained for the choe data are summarized in table  <dig>  for the within-condition analyses the estimates of the proportions of differentially expressed genes are near zero, indicating high specificity of the method.
for the between condition analyses, the number of differentially expressed genes is estimated to be approximately  <dig>  of the genes declared differentially expressed approximately 95% are true positives, demonstrating the methods high positive predictive value.
a similar analysis on the affyu133a data is summarized in figure  <dig>  supplementary figure  <dig>  and table  <dig>  for this data set the histograms of the p values obtained for the replicate array comparisons and the between experiment comparisons are more similar.
with the very few spike-in genes, and thus little difference between the within and between experiment analyses in this data set, this is expected.
also note that the central component is tighter for the affyu133a data set than for the choe data set, reflecting the lower level of noise.
focusing on the tails of the histograms for the analysis of arrays from different experiments , there is a clear deviation from the shape expected under the null, in terms of an excess of p values near  <dig> as well as  <dig>  indicating the presence of both over- and under-expressed genes.
the estimated numbers of differentially expressed genes in the between experiment analyses are around  <dig>  and only a couple for the within experiment analyses.
of the genes declared differentially expressed in the between condition comparisons approximately half are true positives.
thus, in spite of the relatively low fold changes of  <dig> in this data set, the method retains good sensitivity and high specificity.
comparison to other methods
few methods are available for performing differential expression analysis from genechip arrays in the absence of replicates.
here we compare the performance of the bgx based method to those of three other methods that may be applied when only a single replicate is available: the ebarrays method of kendziorski et al.
<cit>  and the affymetrix mas comparison calls  <cit> .
the mas comparison calls are based on the wilcoxon signed rank test applied, for each gene, to the sets of pm-mm values on the two arrays to be compared.
the wilcoxon signed rank test is available in r. we rank the genes on their p-values and use as cut-off the recommended value of  <dig>  to declare the genes as differentially expressed .
the ebarrays method  implements the empirical bayes gamma-gamma or lognormal-normal mixture models, originally developed for two-colour cdna arrays, but equally applicable to analysis of differential expression between genechip arrays.
we use the lognormal-normal model, as generally recommended for genechip data, and apply the method to both the rma expression values  and the mas <dig> values.
the method estimates the proportion of non differentially expressed genes, p <dig>  and for each gene, their posterior probability of belonging to the non-null component.
we rank the genes following decreasing values of this probability and compute the rank cut-off by multiplying p <dig> by the total number of genes analysed.
as an additional comparison, we show results for the efron  method applied to the standardised bgx differences zg = mean/sd, g =  <dig> ...,g.
the method fits a spline to histogram counts of the zg values, using poisson regression.
the null distribution is then empirically estimated from the observed distribution by a first and second moment fitting of a normal distribution to the central component, and used to produce an estimate of p <dig>  the genes are ranked on their absolute zg values and a rank cut-off is estimated as for the ebarrays method above.
estimated numbers of differentially expressed genes obtained with the various methods are given in table  <dig>  along with numbers and proportions of true positives among those declared differentially expressed.
the specificity and positive predictive value is highest for the bgx based method.
among the other methods, the ebarrays method with rma expression measures performs best.
at the opposite extreme is the ebarrays method with mas <dig> expression values: this combination results in very poor specificity and lowest positive predictive value of the methods tried.
the poor performance is likely to be caused by the high variability and strong mean variance relationship of the mas <dig> measure, which in combination with the assumption of a common variance for the genes and constant coefficient of variation of the ebarrays model is detrimental.
in contrast, the rma method has low, stable, variance across the genes, and thus conforms well to the assumptions of the ebarrays method.
efron's method applied to the standardised bgx differences has relatively high positive predictive value .
however, as seen in the large number of wrongly declared differentially expressed genes in the within condition comparisons, this is achieved at the expense of high numbers of false positives.
this shows that the approximation of the null component of the histogram of standardised bgx differences by a gaussian distribution in the efron method is not entirely appropriate and that a more flexible, non-parametric fit is indeed needed.
finally, the performance of the wilcoxon signed rank test is relatively poor on both positive predictive value and specificity.
a biological example
the applications of the method on the spike-in data set analyses illustrate the performance of the method for detecting and estimating the proportion of genes that have been spiked into a common pool of rna at different concentrations.
this is a somewhat stylized situation and, in practice, what we are interested in is the performance of the proposed method in analyses of realistic biological data sets.
the virtue of the spike-in data sets is that we know exactly which genes are differentially expressed.
the problem with testing the method on a biological data set is that there are no biological data sets which possess this feature.
in the absence of such a data set we consider instead a biological data set for which we will regard part of the truth as known: the lin data set  <cit> .
lin et al.
conducted a time course experiment on mice, with the aim of identifying hair cycle-associated genes.
in mice, the first two hair growth cycles are synchronous, but after these, hair growth proceeds asynchronously.
in the lin et al.
study, the expression of genes in mouse back skin was measured at  <dig> time points taken during the first hair-growth cycle and at three time points after the second hair-growth cycle.
the authors hypothesized that genes related to hair-growth cycle should exhibit increased replicate variance from the synchronous to the asynchronous phase, and identified  <dig> genes for which the replicate variance was significantly increased.
to validate their findings they went through a literature search and compiled a list of  <dig> genes that have been shown to be hair cycle-dependent by using other methods .
of the compiled list of  <dig> genes,  <dig> were among the  <dig> identified by lin et al's increased replicate variance based method on the genechip data, thus verifying that hair-cycle associated genes were found by their method.
the time-course profiles of the  <dig> genes were subsequently clustered, the clusters studied and found to relate to distinct genetic pathways.
we will use the literature-based compiled list of  <dig> hair cycle-associated genes as a starting point for making a list of what we will assume to be a subset of the list of "truly differentially expressed genes" between two select time points in the lin data set.
note that a hair-cycle associated gene need not be differentially expressed at any pair of the time points studied.
an examination of the time course cluster profiles indicates that a sizeable subset of the hair-cycle associated genes may have different expression levels at time points 14d and 23d , and we choose to consider arrays from these two time points.
of the  <dig> genes in the literature compiled list, only the  <dig>  that were included in the list of the  <dig> genes were clustered, and we thus only have information on the expression levels at the two considered time points for these  <dig> genes.
the cluster membership of the  <dig> genes suggest that  <dig> belong to clusters with similar levels of expression at time points 14d and 23d.
thus, we are left with a list of  <dig> genes, from the initial list of  <dig>  that:  have been shown to be hair-cycle associated with other methods such as rt-pcr, and that,  we would expect to have different levels of expression at time points 14d and 23d .
importantly, there may well be many other genes on the arrays, apart from those associated to the hair growth cycle associated genes, that are differentially expressed at the two time points considered.
thus our assumption is that the list of  <dig> is a subset of the list  of "truly differentially expressed genes at the two time points".
we take as a criterion for success the ability of the method to detect as differentially expressed these  <dig> genes – allowing that more genes than  <dig> may well be differentially expressed.
CONCLUSIONS
we described a method for performing differential expression analysis from genechip arrays which does not require replicates.
in the method the posterior distributions of the expression differences are used to obtain ranked gene lists and to estimate the proportion of differentially expressed genes.
we investigated the performance of the method for analysis of one single array against another and found the method to perform well on controlled data sets as well as in a biological setting.
the method is not limited to analysis of one single array against another, but is indeed applicable to analyses of any one number of replicates against any other number.
the method for estimating the proportion of differentially expressed genes relies on the empirical estimation of the null distribution of the p values.
for this, we modified the method of efron  <cit>  in two ways.
efron's method operates on standardised differences and fits a normal to the central part of the observed distribution of these values.
we found that the normal assumption for the null component of the standardised bgx differences did not produce convincing results.
instead, we chose to use the natural quantity: p that compares the expression difference to its value under the null, and to exploit its expected shape under the null.
the values of p are bounded by  <dig> and  <dig>  and there is no reason to expect that the symmetric smoothly central component of their distribution should follow any particular parametric shape.
fitting a mixture of normal distributions with a variable number of components as in richardson and green  <cit>  results in the central  component being fitted by a number of normals with means near  <dig> , and different variances.
this makes characterising which components of the mixture make up the null difficult.
in the procedure that we describe we have opted instead for simplicity and attempting to stay close to the empirical distribution by using a spline fit, rather than fitting a particular parametric form.
a different approach to differential expression analysis in the bgx framework would be add an additional hierarchical level to our bgx model and to define a mixture prior for the expression levels rather than a flat prior.
this approach has be taken by a number of authors, for example kendiorski et al.
<cit>  and gottardo et al.
<cit> , working from  expression level data, rather than probe level intensities.
such an extended bgx model would be statistically and theoretically appealing, but also computationally demanding, and will be explored in future work.

 affymetrix genechips are one of the most widely used commercially available oligonucleotide arrays. they have gained widespread popularity for a number of reasons, among which are their high degree of standardization of the production process and their ability to interrogate tens of thousands of genes simultaneously. they differ from many other array types in that they are one color arrays and because each gene is represented by a probe set of pairs of perfect match and mismatch probes. each pm is chosen to match a particular base pair stretch of the sequence encoding the gene, and the full set of pms is chosen with the aim of uniquely identifying the gene. the accompanying mms are identical to their pm counterparts except for a complementary base substitution at the middle nucleotide. they are intended to be used to correct for non-specific hybridization. the full set of pms and mms for a gene represent the basis for the estimation of the level of expression of the gene. most genechip based expression studies are carried out with few replicates. there are two major factors behind this: the considerable cost of the genechip arrays and limitations on the amount of available rna. as microarray expression studies are prone to experimental imperfections, failed hybridizations are often encountered resulting in further reduction in the number of replicates available for analysis. if the number of replicates falls below three, most available analysis tools become unsuitable or unapplicable because they rely on the estimation of variances which is difficult in such circumstances. thus, the development of methods for analyzing experiments with few or no replicates is of high importance. bayesian gene expression, bgx, is an integrated approach to the analysis of affymetrix genechip arrays. it relies on the formulation of a bayesian hierarchical model for estimating expression levels from probe level genechip data. in the bgx approach background correction, gene expression level estimation and differential expression is performed in an integrated analysis, allowing all uncertainties to be dealt with simultaneously in a coherent statistical framework. posterior distributions of gene and condition specific bgx expression levels are obtained from a simultaneous consideration of the probe pair intensities in the available probe sets representing the gene. samples from the posterior distributions are generated using markov chain monte carlo methods. if replicate arrays are available the information in their probe set intensities will be considered jointly in the estimation of expression levels. replicate arrays, however, are not essential for obtaining the posterior distributions of expression for the genes – these will be obtained from the collection of intensities on the array even when only a single array is available. samples from the posterior distribution of the differences in expression levels present a basis for inference on differential expression. in this paper we develop a method for performing differential expression studies from genechip experiments without replicates. the procedure exploits the posterior distributions of differences in expression, obtained from a bgx analysis of genechip arrays, for creating ranked gene lists. in order to define suitable cut-offs for the list, an estimate of the proportion of differentially expressed genes is obtained by empirically estimating the null distribution of a relevant statistic using an approach similar to that of efron. the performance of the method is tested on publicly available spike-in data sets and compared to those of other methods, and further evaluated in a biological study.  the bgx model and methodology most methods for analyzing genechip arrays adopt a stepwise procedure for obtaining a point estimate of expression for each gene on each array. the steps in the procedures typically consist of background correction and normalization followed by summarization or the fitting of a linear model, often performed on the log-scale background corrected intensities. having obtained a point estimate of expression for each gene on each array, studies of differential expression between pairs of conditions are carried out by comparing the collections of point estimates under the conditions using a t-type statistic such as in sam, cyber-t or limma. the bgx method differs from these stepwise point estimate approaches in that all steps in the analysis are dealt with simultaneously, gene and condition specific expression levels are estimated from a joint consideration of the available probe set intensities and the outcomes are posterior distributions of expression rather than point estimates. thus, uncertainties associated with each of the steps are taken into account at all levels of analysis, and the joint uncertainty on the expression level for a gene under a condition is reflected in the shape of the posterior distribution obtained for the level of expression for that gene under that condition. explicitly, with g, j, c and r denoting gene, probe, condition and replicate, respectively, let sgjcr and hgjcr denote gene, probe, condition and replicate specific and non-specific binding of rna, and let φ ∈ be a fraction. to further allow for additive array specific errors, e.g. accommodating mms bigger than pms, the bgx model hypothesizes: pmgjcr ~ n mmgjcr ~ n.     information on the level of expression of gene g under condition c is represented by the set of signal parameters representing the gene under this condition: sgjcr, j =...,j, r =..., rc. we assume that these, shifted and logged, come from a gene and condition specific truncated normal distribution. the non-specific hybridization parameters hgjcr reflect characteristics specific to the sample hybridized, leading us to assume array specific truncated normal distributions for these, shifted and logged. thus, log ~ tn, log ~ tn.     we will here refer to the μgc parameters as the bgx expression indices or levels. we assume exchangeability of the gene and condition specific variance parameters, log() ~ n,     with ac and fixed at values obtained by an empirical bayes like approach, thus stabilizing the variance estimation. in all of the above, the distributions are conditional on variables on the right hand side, and independent for all suffices. the model is fully specified by declaring the following, generally weakly informative priors, independent for all suffices: μgc ~ u, φ ~, λgc ~ n), ~ Γ and ~ Γ). for a more in-depth discussion of the model we refer to. the bgx model relies on mcmc methods for obtaining samples from the posterior distributions of the parameters. the shapes of the posterior distributions of the bgx gene expression indices, μgc, are determined by the probe response patterns. thus, a highly consistent probe set response leads to a tight posterior distribution of expression, and a less consistent pattern will result in a flatter, possibly multi-modal, posterior distribution. examples of posterior expression index distributions, μg, and μg, are given in figure. the corresponding kernel density plots for the differences in expression indices, dg = μg, - μg, are given in figure. the uncertainties of the expression indices are reflected in the shape of these distributions. for gene the multi-modality of the posterior distribution of the expression index under condition is reiterated by the multi-modal posterior distribution of the difference in expression. for the other two example genes the posterior distributions of the differences in expression are tight and uni-modal, centered close to zero and around one respectively, indicating similar expression levels for gene and different expression levels for gene under the two conditions considered. addressing differential expression with replicates a popular approach to conducting differential gene expression studies is to rank the genes according to their degree of evidence for differential expression, and to estimate false discovery rates for different cut-offs on the ranked gene list. this allows the experimenter to obtain a prioritized list of genes to pursue in follow-up studies, with a guidance as to how many genes on the list are expected to be false positives. such approaches are taken in the implementations of the sam, limma and cyber-t methods. each of the methods calculate a different modified t-statistic, the modification relating to the standard deviation or variance calculation in the denominator, and genes are ranked on the resulting p-values. in sam a false discovery rate is estimated based on permuting the original data to get the distribution of t-statistics under the null-hypothesis of no differential expression. limma is implemented with the benjamini and hochberg method for estimating fdr and calculation of adjusted p-values. cyber-t adopts the method of allison et al. for fitting a mixture of beta-distributions distribution) to the observed p-values, and reports estimated true and false positives along with the posterior probability of differential expression. thus, all methods make use of point estimates of expression and depend upon replicates being available for estimation of the variance in the modified t-statistics. addressing differential expression without replicates without replicates the above methods for analysis of differential expression are unapplicable and alternative methods are needed. in this section we describe how the bgx model and methodology may be exploited to obtain such a procedure. we use features of the samples from the posterior distributions collected in the mcmc sampling to produce ranked gene lists. the ranking takes into account the estimated difference in expression level as well as the associated uncertainty. we then consider the set of posterior probabilities of expression differences being smaller than zero, p, g =..., g. by comparing their observed distribution to that expected under the null-hypothesis of no differential expression, we obtain an estimate of the number of differentially expressed genes, gd. this allows us to choose the cut-off of the ranked gene list in an informed manner. the procedures for ranking the genes and for estimating gd in the bgx framework are described in two separate subsections below. the final part of this section contains a comparison of the performance of the bgx based method for performing differential expression analysis from genechips without replicates to those of other available methods: the ebarrays method of kendziorski et al., the wilcoxon signed rank test for comparison calls of the microarraysuite software, and the efron method using the standardized bgx differences as z-statistics. ranking genes using bgx in the bgx framework, the samples from the posterior distributions of the differences dg = μg, - μg, g =...,g, represent a natural base for inference on differential expression between conditions and these are available irrespectively of the number of replicates for each condition used in the analysis. there are numerous ways in which these posterior distributions can be exploited with the aim of addressing expression differences. here we study two types of rankings reflecting the potential of the genes as promising candidates for differential expression: ranking on the'standardized bgx differences', zg = mean/sd, where the mean and standard deviation are computed from the posterior sample of dg values, and ranking on the highest percentile, α*, for which the α-percent credibility interval for the difference dg does not cover zero. note that both rankings use the levels of differential expression as well as the uncertainty of these in the ranking. without replicates point estimate based methods clearly do not have this ability. to illustrate ranking, consider the posterior distributions of expression index differences in figure. all sampled values from the posterior distribution of d are above zero and α* for gene is indistinguishable from 100%. for gene only very tight credibility intervals exclude zero, and α* for this gene is small. for gene α* is around 75%. rankings of type and differ in their emphasis of the posterior distribution characteristics: type summarizes the full distribution by a traditional t-type statistic and we would expect it to perform well when ranking gaussian-shape posterior dg distributions and possibly less well in the presence of asymmetric or multi-modal distributions. type rankings use the tails of the posterior dg distributions and should deal equally well with symmetric and asymmetric distributions. however, we rely on a finite sample from the mcmc procedure for approximating the distributions, and the estimation of the tails of the distributions may be fragile. we examined the performance of the two ranking procedures for the bgx method on the results from the nine analyses of pairs of arrays from the choe data set, that consist of a condition c array and a condition s array. on each pair of arrays we ran the bgx model with conditions and replicate per condition, and obtained ranked gene lists using each of the rankings and. for the choe data it is known exactly which genes are differentially expressed. this allows the exact numbers of true and false positives and negatives to be calculated, for all possible cut-offs in ranks of the ranked gene lists. the two rankings resulted in almost identical counts, and there is no indication from this analysis that either is to be preferred. roc curves summarizing the counts obtained with ranking are shown in figure the performance of the bgx-based rankings from the single replicate comparisons in the choe data set is remarkable: a quarter of the truly differentially expressed genes are included in gene lists with realized false discovery rates of. by extending the gene list to genes the proportion of truly differentially expressed genes detected is increased to 50% and the realized false discovery rate thus about 6%. gene lists of lengths include 70% of the truly differentially expressed genes and have observed false discovery rates of approximately 30%. furthermore, the curves for the nine different analyses of pairs of arrays are highly similar, indicating a stable performance. for comparison, roc curves obtained for the same pairwise analyses using mas and rma are also given in figure for these analyses, genes were ranked on the absolute values of the differences in expression measures between conditions. with just one replicate array per condition, uncertainty of the estimates of expression is not accounted for by these methods, and they both do less well than the bgx-based method. due to the high variability of mas this method performs particularly poorly in the one-versus-one array rankings. to put the above performances into perspective the roc curve obtained for the same method of ranking but using results from a bgx analysis that uses all the available arrays are also shown in figure for this analysis the number of truly differentially expressed genes included in the lists are increased by and 5% respectively to and 75% for the levels of realized false discovery rates of, and. thus, as expected, gene ranking is improved when all replicates are used. however, the proportion of the information contained in the full data set, that may be extracted from the single replicate analyses array is considerable. we note that the performance of the bgx multiple array model on the full choe data set is among the methods found to perform best. results from a similar analysis on the much less extensive affyu133a spike-in data set, are shown in figure as only of the genes represented in this data set are differentially expressed, we plot absolute rather than fractional values of true and false positives. for this data set the ranking of the top genes produced by the three methods are similar. for gene lists longer than bgx and rma outperform mas and rma performs better than bgx for lengths above thus, the relative performances of the methods differ for the choe and the affyu133a data sets. this is most likely due to the different levels of noise in the two data sets. with only spike-in genes in the affyu133a data set arrays from two different conditions are almost like technical replicates. in the choe data set, the genes spiked-in at varying concentrations at small to moderate fold changes result in a noticeable, and more biologically realistic, level of noise between the two conditions compared. thus, for the choe data set, noise has an impact and accounting for this, as is done with the bgx based ranking, is essential, whereas for the affyu133a data set, the impact of noise is negligible and the importance of taking this into account is outweighed by that of reproducibility of point measures for technically replicate arrays. estimation of the proportion of differentially expressed genes having obtained a ranked list of genes the next question is whether we can choose a suitable cut-off. depending on the downstream goal, we may wish to arrive at a list that has a good chance of containing most or all of the "interesting" genes, accepting that false positives will also be included, or we may prefer to end up with a list of genes most of which would appear to be promising candidates for differential expression, expecting very few false positives to be included. to guide the choice of cut-off, it is useful to obtain an estimate of the proportion of differentially expressed genes. to do this, we estimate empirically the distribution of a relevant statistic under the null hypothesis, which in turn allows for the quantification of the proportion of non-null behaving genes, following the idea of efron. we consider, for each gene g, the posterior distribution of the difference in expression, dg = μg, - μg, obtained from the bgx analysis. under the null-hypothesis of no differential expression this posterior distribution should be centered on zero. rephrasing in terms of the posterior p probability, under the null we expect this to be. considering the posterior p probabilities for the full set of genes analyzed, under the null hypothesis of no differential expression there should be decreasingly less support for values away from towards and thus, with no differentially expressed genes a histogram of the p probabilities should be uni-modal with mode of approximately, and have smoothly decreasing tails. the width of the central component of the histogram will depend on the posterior distributions of the dg's: the less clean the information on the expression levels, the less tight the dg distributions, and thus the flatter the histogram of p values. values but could of course equally well have phrased them in terms of the p values). deviations in a histogram of p values from the expected shape under the null hypothesis of no differential expression indicate the presence of differentially expressed genes: excess of p values near zero and one indicate over-expressed and under-expressed genes in condition relative to condition respectively. to quantify the number of non-null genes, we adopt an approach similar to that of efron. we fit a polynomial, f, to the histogram counts by poisson regression and use the central part of the histogram to estimate the null component, while excess area in the tails will represent differentially expressed genes. to be precise, we use the following procedure: we identify the inner global max and left and right most minima on the fitted curve f, and refer to the histogram categories of these points as cmax, and respectively. we use the central part of the histogram to obtain the distribution of the p values under the null: we fit two new curves, and, to the parts of the histogram left and right of category cmax, respectively, using the same fitting procedure as for the curve f but fixing the histogram counts of the outermost categories c and ck to zero, and giving zero weight to categories c..., and,..,ck- respectively. an estimate of the number of genes under the null is obtained by summing the fitted values of the empirical null distribution f denoting the total number of genes by g, the estimated proportion of differentially expressed genes is estimated numbers of over- and under-expressed genes, and, are obtained by quantifying the excess genes in either tail of the histogram relative to those expected under the null. denoting the histogram count in category ci by h we set we applied the above method to the bgx analyses of the pairs ofarrays from the choe and affyu133a data sets. an example of a histogram of p values from a within condition comparison of two arrays from the choe data set is given in figure upper panel. the full set of histograms for within condition analyses of pairs of choe arrays are summarized in figure left, in terms of plots of the curves, f, fitted to the histogram counts by poisson regression. for within condition comparisons the arrays are replicates, so there should be no differentially expressed genes, and the plots indeed exhibit the shape expected under the null hypothesis: they have a central mode near and smoothly decreasing tails. they thus confirm our expectations and indicate that the method works well under the null. figure lower panel, and figure right panel, display the equivalent plots for pairwise between condition analyses of the choe data set arrays. the histograms for these analyses exhibit a clear deviation from the shape expected under the null in terms of an excess of small p values near zero indicating the presence of over-expressed genes. the right-hand tails of the histograms decrease smoothly and there is no indication of under-expressed genes. this is exactly the pattern that should emerge for the choe between condition comparisons: all differentially spiked-in genes have higher concentration under condition s than under condition c. the estimates of the numbers of differentially expressed genes obtained for the choe data are summarized in table for the within-condition analyses the estimates of the proportions of differentially expressed genes are near zero, indicating high specificity of the method. for the between condition analyses, the number of differentially expressed genes is estimated to be approximately of the genes declared differentially expressed approximately 95% are true positives, demonstrating the methods high positive predictive value. a similar analysis on the affyu133a data is summarized in figure supplementary figure and table for this data set the histograms of the p values obtained for the replicate array comparisons and the between experiment comparisons are more similar. with the very few spike-in genes, and thus little difference between the within and between experiment analyses in this data set, this is expected. also note that the central component is tighter for the affyu133a data set than for the choe data set, reflecting the lower level of noise. focusing on the tails of the histograms for the analysis of arrays from different experiments, there is a clear deviation from the shape expected under the null, in terms of an excess of p values near as well as indicating the presence of both over- and under-expressed genes. the estimated numbers of differentially expressed genes in the between experiment analyses are around and only a couple for the within experiment analyses. of the genes declared differentially expressed in the between condition comparisons approximately half are true positives. thus, in spite of the relatively low fold changes of in this data set, the method retains good sensitivity and high specificity. comparison to other methods few methods are available for performing differential expression analysis from genechip arrays in the absence of replicates. here we compare the performance of the bgx based method to those of three other methods that may be applied when only a single replicate is available: the ebarrays method of kendziorski et al. and the affymetrix mas comparison calls. the mas comparison calls are based on the wilcoxon signed rank test applied, for each gene, to the sets of pm-mm values on the two arrays to be compared. the wilcoxon signed rank test is available in r. we rank the genes on their p-values and use as cut-off the recommended value of to declare the genes as differentially expressed. the ebarrays method implements the empirical bayes gamma-gamma or lognormal-normal mixture models, originally developed for two-colour cdna arrays, but equally applicable to analysis of differential expression between genechip arrays. we use the lognormal-normal model, as generally recommended for genechip data, and apply the method to both the rma expression values and the mas values. the method estimates the proportion of non differentially expressed genes, p and for each gene, their posterior probability of belonging to the non-null component. we rank the genes following decreasing values of this probability and compute the rank cut-off by multiplying p by the total number of genes analysed. as an additional comparison, we show results for the efron method applied to the standardised bgx differences zg = mean/sd, g =...,g. the method fits a spline to histogram counts of the zg values, using poisson regression. the null distribution is then empirically estimated from the observed distribution by a first and second moment fitting of a normal distribution to the central component, and used to produce an estimate of p the genes are ranked on their absolute zg values and a rank cut-off is estimated as for the ebarrays method above. estimated numbers of differentially expressed genes obtained with the various methods are given in table along with numbers and proportions of true positives among those declared differentially expressed. the specificity and positive predictive value is highest for the bgx based method. among the other methods, the ebarrays method with rma expression measures performs best. at the opposite extreme is the ebarrays method with mas expression values: this combination results in very poor specificity and lowest positive predictive value of the methods tried. the poor performance is likely to be caused by the high variability and strong mean variance relationship of the mas measure, which in combination with the assumption of a common variance for the genes and constant coefficient of variation of the ebarrays model is detrimental. in contrast, the rma method has low, stable, variance across the genes, and thus conforms well to the assumptions of the ebarrays method. efron's method applied to the standardised bgx differences has relatively high positive predictive value. however, as seen in the large number of wrongly declared differentially expressed genes in the within condition comparisons, this is achieved at the expense of high numbers of false positives. this shows that the approximation of the null component of the histogram of standardised bgx differences by a gaussian distribution in the efron method is not entirely appropriate and that a more flexible, non-parametric fit is indeed needed. finally, the performance of the wilcoxon signed rank test is relatively poor on both positive predictive value and specificity. a biological example the applications of the method on the spike-in data set analyses illustrate the performance of the method for detecting and estimating the proportion of genes that have been spiked into a common pool of rna at different concentrations. this is a somewhat stylized situation and, in practice, what we are interested in is the performance of the proposed method in analyses of realistic biological data sets. the virtue of the spike-in data sets is that we know exactly which genes are differentially expressed. the problem with testing the method on a biological data set is that there are no biological data sets which possess this feature. in the absence of such a data set we consider instead a biological data set for which we will regard part of the truth as known: the lin data set. lin et al. conducted a time course experiment on mice, with the aim of identifying hair cycle-associated genes. in mice, the first two hair growth cycles are synchronous, but after these, hair growth proceeds asynchronously. in the lin et al. study, the expression of genes in mouse back skin was measured at time points taken during the first hair-growth cycle and at three time points after the second hair-growth cycle. the authors hypothesized that genes related to hair-growth cycle should exhibit increased replicate variance from the synchronous to the asynchronous phase, and identified genes for which the replicate variance was significantly increased. to validate their findings they went through a literature search and compiled a list of genes that have been shown to be hair cycle-dependent by using other methods. of the compiled list of genes, were among the identified by lin et al's increased replicate variance based method on the genechip data, thus verifying that hair-cycle associated genes were found by their method. the time-course profiles of the genes were subsequently clustered, the clusters studied and found to relate to distinct genetic pathways. we will use the literature-based compiled list of hair cycle-associated genes as a starting point for making a list of what we will assume to be a subset of the list of "truly differentially expressed genes" between two select time points in the lin data set. note that a hair-cycle associated gene need not be differentially expressed at any pair of the time points studied. an examination of the time course cluster profiles indicates that a sizeable subset of the hair-cycle associated genes may have different expression levels at time points 14d and 23d, and we choose to consider arrays from these two time points. of the genes in the literature compiled list, only the that were included in the list of the genes were clustered, and we thus only have information on the expression levels at the two considered time points for these genes. the cluster membership of the genes suggest that belong to clusters with similar levels of expression at time points 14d and 23d. thus, we are left with a list of genes, from the initial list of that: have been shown to be hair-cycle associated with other methods such as rt-pcr, and that, we would expect to have different levels of expression at time points 14d and 23d. importantly, there may well be many other genes on the arrays, apart from those associated to the hair growth cycle associated genes, that are differentially expressed at the two time points considered. thus our assumption is that the list of is a subset of the list of "truly differentially expressed genes at the two time points". we take as a criterion for success the ability of the method to detect as differentially expressed these genes – allowing that more genes than may well be differentially expressed.  we described a method for performing differential expression analysis from genechip arrays which does not require replicates. in the method the posterior distributions of the expression differences are used to obtain ranked gene lists and to estimate the proportion of differentially expressed genes. we investigated the performance of the method for analysis of one single array against another and found the method to perform well on controlled data sets as well as in a biological setting. the method is not limited to analysis of one single array against another, but is indeed applicable to analyses of any one number of replicates against any other number. the method for estimating the proportion of differentially expressed genes relies on the empirical estimation of the null distribution of the p values. for this, we modified the method of efron in two ways. efron's method operates on standardised differences and fits a normal to the central part of the observed distribution of these values. we found that the normal assumption for the null component of the standardised bgx differences did not produce convincing results. instead, we chose to use the natural quantity: p that compares the expression difference to its value under the null, and to exploit its expected shape under the null. the values of p are bounded by and and there is no reason to expect that the symmetric smoothly central component of their distribution should follow any particular parametric shape. fitting a mixture of normal distributions with a variable number of components as in richardson and green results in the central component being fitted by a number of normals with means near, and different variances. this makes characterising which components of the mixture make up the null difficult. in the procedure that we describe we have opted instead for simplicity and attempting to stay close to the empirical distribution by using a spline fit, rather than fitting a particular parametric form. a different approach to differential expression analysis in the bgx framework would be add an additional hierarchical level to our bgx model and to define a mixture prior for the expression levels rather than a flat prior. this approach has be taken by a number of authors, for example kendiorski et al. and gottardo et al., working from expression level data, rather than probe level intensities. such an extended bgx model would be statistically and theoretically appealing, but also computationally demanding, and will be explored in future work. 
we exploit these posterior distributions to create ranked gene lists that take into account the estimated expression difference as well as its associated uncertainty. we estimate the proportion of non-differentially expressed genes empirically, allowing an informed choice of cut-off for the ranked gene list, adapting an approach proposed by efron. the bgx method estimates a posterior distribution of expression for each gene and condition, from a simultaneous consideration of the available probe intensities representing the gene in a condition. importantly, posterior distributions of expression are obtained regardless of the number of replicates available. we assess the performance of the method, and compare it to those of other methods, on publicly available spike-in data sets, as well as in a proper biological setting. studies of differential expression that use affymetrix genechip arrays are often carried out with a limited number of replicates. the procedure relies on a bayesian integrated approach to the analysis of affymetrix genechips. we describe a statistical procedure for performing differential expression analysis without replicates. most existing methods for studying differential expression rely on the availability of replicates and the demand for alternative methods that require few or no replicates is high. the method presented is a powerful tool for extracting information on differential expression from genechip expression studies with limited or no replicates.
 studies of differential expression that use affymetrix genechip arrays are often carried out with a limited number of replicates. reasons for this include financial considerations and limits on the available amount of rna for sample preparation. in addition, failed hybridizations are not uncommon leading to a further reduction in the number of replicates available for analysis. most existing methods for studying differential expression rely on the availability of replicates and the demand for alternative methods that require few or no replicates is high.  we describe a statistical procedure for performing differential expression analysis without replicates. the procedure relies on a bayesian integrated approach to the analysis of affymetrix genechips. the bgx method estimates a posterior distribution of expression for each gene and condition, from a simultaneous consideration of the available probe intensities representing the gene in a condition. importantly, posterior distributions of expression are obtained regardless of the number of replicates available. we exploit these posterior distributions to create ranked gene lists that take into account the estimated expression difference as well as its associated uncertainty. we estimate the proportion of non-differentially expressed genes empirically, allowing an informed choice of cut-off for the ranked gene list, adapting an approach proposed by efron. we assess the performance of the method, and compare it to those of other methods, on publicly available spike-in data sets, as well as in a proper biological setting.  the method presented is a powerful tool for extracting information on differential expression from genechip expression studies with limited or no replicates. 