A class of models for analyzing GeneChip^{®} gene expression analysis array data
 Wenhong Fan^{1},
 Joel I Pritchard^{2},
 James M Olson^{2},
 Najma Khalid^{1} and
 Lue Ping Zhao^{1}Email author
DOI: 10.1186/14712164616
© Fan et al; licensee BioMed Central Ltd. 2005
Received: 28 July 2004
Accepted: 14 February 2005
Published: 14 February 2005
Abstract
Background
Various analytical methods exist that first quantify gene expression and then analyze differentially expressed genes from Affymetrix GeneChip^{®} gene expression analysis array data. These methods differ in the choice of probe measure (quantification of probe hybridization), summarizing multiple probe intensities into a gene expression value, and analysis of differential gene expression. Research papers that describe these methods focus on performance, and how their approaches differ from others. To better understand the common features and differences between various methods, and to evaluate their impact on the results of gene expression analysis, we describe a class of models, referred to as generalized probe models (GPMs), which encompass various currently available methods.
Results
Using an empirical dataset, we compared different formulations of GPMs, and GPMs with three other commonly used methods, i.e. MAS 5.0, dChip, and RMA. The comparison shows that, on a genomewide scale, different methods yield similar results if the same probe measures are chosen.
Conclusion
In this paper we present a general framework, i.e. GPMs, which encompasses various methods. GPMs permit the use of a wide range of probe measures and facilitate appropriate comparison between commonly used methods. We demonstrate that the dissimilar results stem primarily from different choice of probe measures, rather than other factors.
Background
Microarray experiments are routinely conducted to assess associations of experimental factors (or disease outcomes) with gene expression profiles. The Affymetrix GeneChip^{®} gene expression analysis array, one of most commonly used microarray technologies, uses multiple oligonuleotides (25mers) to measure expression abundance of a single gene. Recognizing that nonspecific hybridization could significantly alter the accurate quantification of transcript abundance, Affymetrix designs the array to contain two types of probes. Probes that are perfectly complementary to the target sequence, called Perfect Matches (PM), are intended to measure mainly specific hybridization. A second set of probes identical to PM except for a single nucleotide in the center of the probe sequence (the 13^{th} nucleotide), called Mismatches (MM), are intended to quantify nonspecific hybridization [1]. A PM and its corresponding MM constitutes a probe pair, and multiple probe pairs, i.e. a probe set, are summarized to measure transcript abundance for a particular gene. "Probe measure" is used in this paper to refer to the manner in which probe hybridization is quantified based on a pair of PM and MM intensity values. For example, PMMM is a probe measure, and PM only is another probe measure.
A number of methods have been developed to quantify gene expression abundance from GeneChip^{®} expression analysis array data using different probe measures and summary schemes. Among them, Microarray Suite 5.0 (MAS 5.0) [1], dChip [2] and robust multiplearray average (RMA) [3] are the best known.
Prior to MAS 5.0, the probe measure used in MAS 4.0 was PMMM [4]. The problem arises when a significant proportion of MM values, (~33% in the HuGeneFL array and ~25% in the Human Genome U133A array), is greater than the corresponding PM values, which makes PMMM negative. To resolve this anomaly, in MAS 5.0, Affymetrix computes an "ideal mismatch" (IM) based on missing data theory such that PMIM is always greater than zero [1]. Then, all probe pairs are used to estimate a gene expression value based on Tukey's Biweight algorithm. However, even with the use of IM, the variation among probes could be greater than between samples.
Li and Wong modelled probe level data to generate model based expression index (MBEI) and implemented it in the dChip software [2]. Noting that probe specificity is significant, highly reproducible and predictable, Li and Wong used a hybridization rate parameter to account for the hybridization specificity for a probe. For a probe pair, hybridization rates are different for PM and MM; the former is always greater than the latter, and both are greater than zero. The rate was fixed for the same probe across all the samples. Both PM and MM together or PM only, can be used in the Li and Wong model.
Another approach, RMA, available from Bioconductor [5], summarizes probe intensities into a gene expression measure based on an additive model on the logarithmic scale of a background corrected PM (PM_{rma}) [3]. RMA estimates a common mean nonspecific hybridization background (for an entire chip) from PM using a convolution model and then subtracts this background from PM to generate the PM_{rma}.
The gene expression obtained from either MAS 5.0 or dChip or RMA can then be used to associate the gene expression values with experimental factors using an algorithm of the users' choice. Three main factors affect the analytical results of differential gene expression analysis: the probe measure chosen, the algorithm used to summarize probe level data into gene expression (called summary algorithm in this paper), and the model used to associate gene expression with the experimental factors (called association model). Direct comparisons of the various approaches proposed for analyzing GeneChip^{®} gene expression data are complicated considering these three factors. Generalizing the various algorithms into one framework would facilitate comparisons.
In this paper we propose a class of generalized probe models (GPMs) that includes various analytical approaches for GeneChip^{®} gene expression analysis array data as special cases. Using an empirical dataset, we assess the impact of different processes on the analytical results by comparing different formulations of GPM as well as GPMs with three other methods, MAS 5.0, dChip, and RMA.
Results
We applied GPM to the analyses of data obtained from a study investigating gene response to ATRA (all trans retinoic acid) or drug diluent (ETOH, ethyl alcohol). Briefly, at twentyfour hours after treatment, total RNA was extracted from cells, processed and hybridized to the HuGeneFL GeneChip^{®}. The dataset consists of ten samples from the ATRA treatment group and ten samples from the control group (ETOH treated) in four medulloblastoma cell lines. We are interested in identifying genes that are differentially expressed between the two treatment groups.
We used three different probe measures: PMIM, PM only and PM_{rma}, and compared the performance of different methods using standardized coefficients, defined as the estimated coefficient divided by its standard error. The reason for using this index is that the standardized coefficients, usually known as Zscore test statistics, are independent of scale, and may be used to make statistical inference.
GPM1 (2), GPM2 (3) and GPM3 (4) can be derived from the full GPM model (1) by making different statistical modeling assumptions. GPM1 takes summarized gene expressions and associates them with experimental factors; GPM2 and GPM3 directly associate probe level data with experimental factors without first summarizing gene expressions (see Methods).
Comparison of GPMs with three other commonly used methods
Comparison among the GPMs
Estimated parameters for eight candidate genes from GPMs
Probe set ID  Y00291_at  L13738_at  D79990_at  M13666_at  X02158_rna1_at  X84002_at  L19605_at  M60503_at  

Parameters  β  SE  Z  β  SE  Z  β  SE  Z  β  SE  Z  β  SE  Z  β  SE  Z  β  SE  Z  β  SE  Z 
GPM1 PMIM  3.17  0.48  6.65  0.80  0.18  4.50  3.04  0.68  4.47  1.07  0.25  4.24  1.07  0.25  4.24  0.63  0.15  4.18  0.51  0.12  4.16  2.05  0.50  4.13 
GPM2 PMIM  2.79  0.36  7.74  0.65  0.13  4.89  2.92  0.67  4.3  1.11  0.24  4.6  0.92  0.24  3.86  0.33  0.12  2.63  0.51  0.09  5.53  1.56  0.35  4.44 
GPM3 PMIM  3.05  0.35  8.60  0.71  0.12  5.69  3.10  0.63  4.89  1.31  0.22  5.93  1.40  0.25  5.73  0.25  0.11  2.26  0.52  0.09  6.02  2.25  0.38  5.96 
In summary, we conclude that the GPMs are similar to MAS 5.0, dChip and RMA on a genomewide level when using the same probe measures and that the choice of probe measure may be more important than the summary algorithms to obtain the gene expression or models used to compute the coefficients.
Discussion
In this paper, we have described a general framework that can be used to compare various methods and evaluate their similarities and differences. We found that various methods tend to generate similar results, on a genomewide scale, when the same probe measure is chosen, and probe measure seems to have greater impact on the analytical results than other factors.
In Figure 1, we compared the standardized coefficients estimated with GPM1 using gene expression computed from MAS 5.0, dChip and RMA with their own dictated probe measures. Since we consistently used GPM1 as the association modeling machinery for each analysis, we assessed the combined impact of probe measure and summary algorithm. We found that the results obtained from dChip PM and RMA PM_{rma} were similar to each other, but different from those obtained from MAS 5.0 PMIM. Although dChip PM and RMA PM_{rma} use different summary algorithms, their analytical results are similar due to the PM based probe measures used in both analyses.
In Figure 2, we see again, that on a genomewide scale, results from MAS 5.0, GPM2 and GPM3 are more similar when same probe measure is used than when the probe measures are different, indicating that the probe measure plays a key role in determining the similarity of results from two methods. Our preliminary analyses suggest that the choice of probe measure has bigger impact on the results than summary algorithm and association modeling.
For the three variants of GPMs, we compared the standardized coefficients from GPM2 or GPM3 with those from GPM1 using the gene expression values computed in MAS 5.0, dChip and RMA. From the high R values (and correspondingly, low MSEs) in the six plots shown in Figure 3, we infer that the standardized coefficients obtained from variants of GPMs are similar when they used the same probe measure.
For seven of the eight candidate genes selected by GPM1 using gene expression values generated by MAS 5.0, the genespecific regression coefficients were similar among the MAS 5.0 PMIM, GPM2 PMIM and GPM3 PMIM. This indicates that for these seven genes it makes little difference between using summary measures or modeling directly at the probe level data in GMP2 or GPM3, when the same probe measure is used.
In addition to the three factors we mentioned (i.e. choice of probe measure, summary algorithm and association modeling) that have an impact on analytical results, data preprocessing/normalization could also affect the analytical results. Some researchers combine the probe measure and preprocessing normalization together. Normalization matters the most when the arrays in an experiment are not comparable to each other. In such cases, normalization process could significantly impact on the result. In our case, we normalized the data in the GPMs using a regressionbased approach [6], either at probe level in GPM2 and GPM3, or on gene expression level in GPM1. The expression measures obtained from dChip and RMA were normalized by their own normalization schemes. However, even with the different normalization schemes, probe measure appears to be the primary factor to impact the results in our data set.
A list of selected probe measures
Scenario  Calculation  Annotation 

1. MAS4.0equivalent  Z_{ jik }= (y_{ji 1k} y_{ji 0k})  Direct difference between PM and MM 
2. MAS5.0equivalent  Z_{ jik }= (y_{ji 1k} )  is the Idealized Mismatch (IM) 
3. PM only  Z_{ jik }= y_{ji 1k}  Ignore MM 
4. RMAequivalent  Z_{ jik }= log(y_{ji 1k} )  is the mean background estimated from PM for the k^{th} chip 
5. Log ratio  Z_{ jik }= ln(y_{ji 1k}/ y_{ji 0k})  Difference on the logarithmic scale 
6. Log difference  Z_{ jik }= ln(y_{ji 1k} )  Difference on the logarithmic scale 
7. Log PM  Z_{ jik }= ln(y_{ji 1k})  PM only on the logarithmic scale 
8. BoxCox on PM  Z_{ jik }= (  1) / ω  BoxCox transformation on PM only 
9. BoxCox on PMIM  Z_{ jik }= [(y_{ji 1k} )^{ ω } 1] / ω  BoxCox transformation on PMIM 
To facilitate the evaluation and use of GPMs, we have developed a software program, called ProbePlus that implements our GPMs. This program will be made available to academic researchers through the website http://qge.fhcrc.org/probeplus.
Conclusions
In this paper we describe a general framework to analyze GeneChip^{®} gene expression analysis array data. This framework is flexible to permit comparisons of different methods with respect to the choice of probe measure and analytical models used. We found that different methods yield similar result when probe measures are the same.
Methods
The generalized probe model
A typical probelevel data generated from GeneChip^{®} gene expression analysis array
Sample ID  Probe  PM (1)  1  2  ...  k  ...  K 

Covariate  ID  MM(0)  x _{1}  x _{2}  ...  x _{ k }  ...  x _{ K } 
ORF_{1}  1  1  y _{1111}  y _{1112}  y _{111k}  ...  y _{111K}  
1  0  y _{1101}  y _{1102}  y _{110k}  ...  y _{110K}  
2  1  y _{1211}  y _{1212}  y _{121k}  ...  y _{121K}  
2  0  y _{1201}  y _{1202}  y _{120k}  ...  y _{120K}  
...  
N  1  y _{1N 11}  y _{1N 12}  y _{1N 1k}  ...  y _{1N 1K}  
N  0  y _{1N 01}  y _{1N 02}  y _{1N 0k}  ...  y _{1N 0K}  
...  
ORF_{ j }  1  1  y _{j 111}  y _{j 112}  y _{j 11k}  ...  y _{j 11K}  
1  0  y _{j 101}  y _{j 102}  y _{j 10k}  ...  y _{j 10K}  
2  1  y _{j 211}  y _{j 212}  y _{j 21k}  ...  y _{j 21K}  
2  0  y _{j 201}  y _{j 202}  y _{j 20k}  ...  y _{j 20K}  
...  
i  1  y _{ji 11}  y _{ji 12}  y _{ji 1k}  ...  y _{ji 1K}  
i  0  y _{ji 01}  y _{ji 02}  y _{ji 0k}  ...  y _{ji 0K}  
...  
N  1  y _{jN 11}  y _{jN 12}  y _{jN 1k}  ...  y _{jN 1K}  
N  0  y _{jN 01}  y _{jN 02}  y _{jN 0k}  ...  y _{jN 0K}  
...  
ORF_{ J }  1  1  y _{J 111}  y _{J 112}  y _{J 11k}  ...  y _{J 11K}  
1  0  y _{J 101}  y _{J 102}  y _{J 10k}  ...  y _{J 10K}  
2  1  y _{J 211}  y _{J 212}  y _{J 21k}  ...  y _{J 21K}  
2  0  y _{J 201}  y _{J 202}  y _{J 20k}  ...  y _{J 20K}  
...  
N  1  y _{JN 11}  y _{JN 12}  y _{JN 1k}  ...  y _{JN 1K}  
N  0  y _{JN 01}  y _{JN 02}  y _{JN 0k}  ...  y _{JN 0K} 
In a typical experiment as described above, it is frequently of interest to discover genes that are significantly associated with one or more experimental covariates x_{ k }. For example, consider an experiment to discover genes that are differentially expressed between two groups, x_{ k }takes a binary values: x_{ k }= 0 for the control group and x_{ k }= 1 for the treatment group. To achieve the scientific objective, the analytic procedure is to assess associations of = (Z_{j 1k}, Z_{j 2k},...,Z_{ jNk })' with covariates x_{ k }via the distribution function f (  x_{ k }). In essence, Z_{ jik }are treated as vectors of multivariate correlated outcome variables, and used to identify the probes/genes that are differentially expressed. Recognizing the high dimensionality of multiple probes and multiple genes, we propose to apply a marginal model that uses marginal means to describe relationship of probes/genes with the covariates without the necessity of specifying the full distribution f (  x_{ k }). Our framework directly associates experimental factors with probe intensities and is referred to as the generalized probe model or GPM. We propose the following (1) to describe relationship between and x_{ k },
where (δ_{ k }, λ_{ k }) are chipspecific heterogeneity factors for k th chip [7], τ_{ ji }are gene and probespecific parameters quantifying the mean intensity value for the i th probe of the j th gene, β_{ ji }quantifies the gene and probespecific parameters quantifying the difference between treated and control groups, and v_{ jik }quantifies expression values for individual probe pairs. Lastly, (ξ_{j 1k}, ξ_{j 2k},...,ξ_{ jNk }) represents a vector of gene and probespecific random variations across K independent samples. Since probe pairs are selected to target the j th gene and are spatially arranged by a preselected design to eliminate common artifacts, they may be correlated because of crosshybridizations or spatial dependencies. From the biological perspective, specifying a joint distribution for (ξ_{j 1k}, ξ_{j 2k},...,ξ_{ jNk }) would be difficult, if not impossible. It is thus preferable to leave it unspecified.
The above GPM (1) includes a range of more simplified models based on specific assumptions. First, under the assumption that all probespecific parameters are the same, i.e., τ_{ ji }= τ_{ j }and β_{ ji }= β_{ j }, the general model (1) simplifies to the following model:
and is equivalent to using a summarized gene expression to associate with the experimental factors [7]. For simplicity and comparison with other special models, we refer this model as GPM1.
If one postulates that all probespecific parameters are not the same, but follow an additive probe model, then general model (1) under modeling assumption that β_{ ji }= β_{ j }, with probespecific values (τ_{j 1}, τ_{j 2},...,τ_{ jN }) may be written as
in which estimating β_{ j }is of primary interest. This variation of the general model is referred to as GPM2.
On the other hand, the probe parameters may follow a multiplicative model (in the spirit of Li and Wong's model), then the third model, referred to as GPM3, is derived under the assumption that τ_{ ji }≈ φ_{ ji }τ_{ j }and β_{ ji }≈ φ_{ ji }β_{ j }, and may be written as
where φ_{ ji }denote the multiplicative probespecific effects and can be uniquely determined by constraining the mean to be one.
Estimation and inference
Our estimation procedures do not require any assumptions with respect to the error distribution, since any distributional assumptions, which may be appropriate for some genes, are likely to be violated for other genes. To ensure the robustness of statistical inference, we propose to use generalized estimating equation theory, which has been fully described in a seminal paper by [8]. In the current context, we choose the "working independence" assumption for modeling dependencies between probes [8, 9]., to avoid making any assumptions on dependence structures. The asymptotic variance matrix is estimated with the usual "sandwich" estimator [8]. Diagonal elements in the variance matrix are estimates of marginal variances for all estimated parameters, and are denoted by for the estimated parameters in the model. Both estimates can be used to construct test statistics, such as the ratio of over , known as Waldstatistic. Under the null hypothesis, each statistic has an asymptotic normal distribution when the sample size is sufficiently large, and can therefore be used for making statistical inferences. When the sample size is small, this quantity is treated as a standardized regression coefficient.
List of abbreviations
 ATRA:

All Trans Retinoic Acid
 ETOH:

Ethyl Alcohol
 GPM:

Generalized Probe Model
 IM:

Ideal Mismatch
 MAS:

Microarray Suite
 MM:

Mismatch
 ORF:

Open Reading Frame
 PM:

Perfect Match
 RMA:

Robust Multiple Array average
Declarations
Acknowledgements
We wish to thank Andrew Strand for his helpful discussion during the preparation of this manuscript. This work was supported by grants from National Institutes of Health.
Authors’ Affiliations
References
 Affymetrix: Microarray Suite 5.0 User's Guide. 2001Google Scholar
 Li C, Wong WH: Modelbased analysis of oligonucleotide arrays: expression index computation and outlier detection. Proc Natl Acad Sci U S A. 2001, 98: 316. 10.1073/pnas.011404098.PubMedPubMed CentralView ArticleGoogle 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: 24964. 10.1093/biostatistics/4.2.249.PubMedView ArticleGoogle Scholar
 Affymetrix: Microarray Suite 4.0 User's Guide. 1999Google Scholar
 Bioconductor project. [http://www.bioconductor.org]
 Thomas JG, Olson JM, Tapscott SJ, Zhao LP: An efficient and robust statistical modeling approach to discover differentially expressed genes using genomic expression profiles. Genome Res. 2001, 11: 122736. 10.1101/gr.165101.PubMedPubMed CentralView ArticleGoogle Scholar
 Zhao LP, Prentice R, Breeden L: Statistical modeling of large microarray data sets to identify stimulusresponse profiles. Proc Natl Acad Sci U S A. 2001, 98: 56316. 10.1073/pnas.101013198.PubMedPubMed CentralView ArticleGoogle Scholar
 Liang KY, Zeger SL: Longitudinal data analysis using generalized linear models. Biometrika. 1986, 73: 1322.View ArticleGoogle Scholar
 Prentice RL, Zhao LP: Estimating equations for parameters in means and covariances of multivariate discrete and continuous responses. Biometrics. 1991, 47: 82539.PubMedView ArticleGoogle 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.