General aspects

Microarray assays involve multiple sources of systematic effects (e.g. differences among slides and dye labeling and spatial variation within slides), such that their statistical analysis generally requires a data preprocessing step, usually referred to as normalization, for removing non-biological effects contributing to variation on expression measures. A review and comparison of normalization techniques is provided, e.g., in Jaffrezic et al. (Reference Jaffrezic, De Koning, Boettcher, Bonnet, Buitenhuis, Closset, Déjean, Delmas, Detilleux, Dovc, Duval, Foulley, Hedegaard, Hornshoj, Hulsegge, Janss, Jensen, Jiang, Lavric, Lê Cao, Lund, Malinverni, Marot, Nie, Petzl, Pool, Robert-Granié, SanCristobal, van Schothorst, Schuberth, Sorensen, Stella, Tosser-Klopp, Waddington, Watson, Yang, Zerbe and Seyfert2007).

Normalization approaches can be based either on control genes, or on data-driven selected genes seemingly equally expressed across experimental groups, or most commonly, on all genes represented in the slides. In this latter case, robust statistical techniques are used, with the assumption that the majority of the genes are not differentially expressed. In some cases, e.g., with boutique arrays developed specifically for studies with certain tissue(s) and developmental stage(s), this assumption does not hold and alternative normalization techniques are required (Oshlack et al., Reference Oshlack, Emslie, Corcoran and Smyth2007). Hence, it is apparent that biological knowledge should be used to guide researchers to choose an appropriate data normalization strategy for their experiments.

Other examples of use of prior biological information in preprocessing of microarray data can be found, among others, in Tuikkala et al. (Reference Tuikkala, Elo, Nevalainen and Aittokallio2006) who used information on functional similarities, assessed by semantic correspondences originating from gene ontology (GO) annotations, to improve the selection of relevant genes for missing data imputation, and in Kim et al. (Reference Kim, Doan, Duggal and Bailey-Wilson2007) who investigated the incorporation of genetic relationships among animals for the normalization of expression data in genetical genomics studies.

After this first step of data preprocessing, a higher-level statistical analysis is performed. Depending on the objective(s) of the study such analysis may involve, e.g., clustering, classification techniques or hypothesis testing for detection of differential expression (for a comparison of different approaches, see, e.g., Jaffrezic et al., Reference Jaffrezic, De Koning, Boettcher, Bonnet, Buitenhuis, Closset, Déjean, Delmas, Detilleux, Dovc, Duval, Foulley, Hedegaard, Hornshoj, Hulsegge, Janss, Jensen, Jiang, Lavric, Lê Cao, Lund, Malinverni, Marot, Nie, Petzl, Pool, Robert-Granié, SanCristobal, van Schothorst, Schuberth, Sorensen, Stella, Tosser-Klopp, Waddington, Watson, Yang, Zerbe and Seyfert2007 and Sorensen et al., Reference Sorensen, Bonnet, Buitenhuis, Closset, Déjean, Delmas, Duval, Tosser-Klopp, Hedegaard, Hornshoj, Hulsegge, Jaffrézic, Jensen, Jiang, de Koning, Lê Cao, Nie, Petzl, Pool, Robert-Granié, San Cristobal, van Schothorst, Schuberth, Seyfert, Waddington, Watson, Yang, Zerbe and Lund2007). In this manuscript, we will specifically discuss the comparison of expression profiles across groups or populations using ANOVA models. Such models (including those with random effects) are the most popular in the microarray literature (see, e.g., Jaffrezic et al., Reference Jaffrezic, De Koning, Boettcher, Bonnet, Buitenhuis, Closset, Déjean, Delmas, Detilleux, Dovc, Duval, Foulley, Hedegaard, Hornshoj, Hulsegge, Janss, Jensen, Jiang, Lavric, Lê Cao, Lund, Malinverni, Marot, Nie, Petzl, Pool, Robert-Granié, SanCristobal, van Schothorst, Schuberth, Sorensen, Stella, Tosser-Klopp, Waddington, Watson, Yang, Zerbe and Seyfert2007), due to their flexibility and ease of use as well as the availability of software for their implementation, e.g. MAANOVA and LIMMA, both developed in R language and available at Bioconductor website (http://www.bioconductor.org/).

The ANOVA approach for analyzing microarray data was first presented by Kerr and Churchill (Reference Kerr and Churchill2001), with a model described as:

(1)

where y _adgvr represents the expression intensities on the log scale; μ is an overall constant; A_a, D_d and (AD)_ad are ‘global factors’ that account for variation between arrays and dyes; G_g is the gene effect; (AG)_ag represents array by gene interactions; (DG)_dg is gene-specific dye effects; (VG)_vg is the quantity of interest, related to differential expression of gene g specifically attributable to variety (treatment) v; and ε _adgvr is random residual term with variance σ _ε². Wolfinger et al. (Reference Wolfinger, Gibson, Wolfinger, Bennett, Hamadeh, Bushel, Afshari and Paules2001) extended this approach to allow for random effects as well for gene-specific variance components. For computational convenience, a two-step methodology was proposed. In the first stage, a model including only global factors is fitted to all data simultaneously, which can be written as:

(2)

where the terms are defined as in model

. The second stage uses the estimated residuals from model

, and fits a series of gene specific models, expressed as:

(3)

where μ _g, A_ag, D_dg and V_vg are gene-specific constant, array effects, dye effects and treatment effects, respectively.

The mixed effects modeling approach of Wolfinger et al. (Reference Wolfinger, Gibson, Wolfinger, Bennett, Hamadeh, Bushel, Afshari and Paules2001) relaxes the assumption of a common variance for all genes, as in model (1), but at the expense of performing independent analyses for each gene, one at a time. This means that genes are treated as independent entities, ignoring the natural covariance among their expression levels due to coregulation processes. The main drawback of this approach is that, because sample sizes for each gene are usually relatively small, the resulting statistical tests have low power and fold change estimates have low precision. Most of the strategies available in the literature to overcome this situation refer to some sort of shrinkage estimation of variance components, to borrow information across genes (e.g. Smyth, Reference Smyth2004 and Cui et al., Reference Cui, Hwang, Qiu, Blades and Churchill2005). However, those approaches are still based on the assumption of independence among genes, and they do not take advantage of prior genomic or biological knowledge. Exogenous information is generally used only after a list of genes has been selected, e.g., to study gene set enrichment (Huang et al., Reference Huang, Sherman and Lempicki2009).

Modeling covariance between gene expression

Some recent publications have suggested the integration of prior biological knowledge into the statistical analysis of microarray data. For example, Imoto et al. (Reference Imoto, Higuchi, Goto, Kuhara, Miyano and Tashiro2003) and Werhli and Husmeier (Reference Werhli and Husmeier2007) exploit multiple sources of prior knowledge for gene regulatory network reconstruction; and Huang and Pan (Reference Huang and Pan2006), Pan (Reference Pan2006a) and Tseng (Reference Tseng2007) propose biologically informed gene clustering techniques. However, not much research has been published regarding the incorporation of prior knowledge into the analysis of differential expression, one exception being Pan (Reference Pan2006b), who uses gene functional annotation for improving detection of differentially expressed genes. Here, we discuss some strategies for integrating prior biological information into the statistical analysis of differential expression. The main idea is to build parsimonious multivariate ANOVA models, as an alternative to the shrinkage estimation approaches usually found in the literature.

We start by combining the gene-specific models into a single multigene model (as discussed by Hoeschele and Li, Reference Hoeschele and Li2005), which can be expressed in matrix notation as:

where y is a vector with expression intensities for all genes and samples; β and u are vectors with fixed (such as treatment and dye) and random (such as slide and animal) effects, respectively; X and Z are incidence matrices; and ε is a vector of residual terms, assumed to be independent from each other. For modeling the dependence across genes, a multivariate version of the mixed model described above could be considered, with an unstructured covariance matrix for modeling the residuals. In this case, the variance of ε could be represented as Var (ε) = I_n⊗Σ, where

is a matrix describing the variances and covariances between genes, k is the number of genes and I_n is an identity matrix of order n (n is the number of samples or subjects in the experiment).

The covariance matrix Σ, however, has k(k + 2)/2 distinct elements that need to be estimated (k variances and k(k − 1)/2 covariances), which makes fitting such a model unfeasible given the small number of samples relative to the number of genes, i.e. n ≪ k. For example, with an experiment involving k = 10 000 genes, there would be over 50 million (co)variance parameters to be estimated. Nonetheless, currently microarray experiments in livestock do not involve more than a few hundreds of slides.

To overcome this problem, biological information could be used, e.g., to split the genes into subsets of potentially coexpressed genes such that multivariate models of much lower dimensions could be fitted within each set of genes. Biologically meaningful gene groupings, such as regulatory networks and pathways, can be assessed using freely available databases such as the GO Consortium (GO – http://www.geneontology.org/), the Kyoto Encyclopedia of Genes and Genomes (KEGG – http://www.genome.jp/kegg/) among others. In this context, the assumption is that genes in different groups are independent but that genes within each group can have any pattern of coexpression. Under these circumstances, the new covariance matrix Σ is block diagonal, and can be described as:

where b is the number of gene subsets, and

The number of non-zero parameters in Σ is now which is much smaller than the original parameters, where . Moreover, the analysis under such settings can be implemented independently for each gene set, making model dimensions even smaller and computations less cumbersome. Nonetheless, small samples sizes can still hinder the gene set models to be fitted, especially for larger groups of genes. In this case, parsimonious structures for each covariance matrix B_i can be built using biological information, e.g., by describing covariances between genes as a function of their number of shared elements, their positions in the genome, as well as information on biological processes (Vazquez et al., Reference Vazquez, Rosa, de los Campos, Weigel and Gianola2009).

A standard way of summarizing current knowledge about the biological process is via a regulatory network. Wei and Pan (Reference Wei and Pan2008), e.g., used biological information summarized in a network with a model in which the joint prior probability of gene effects is multivariate normal, and the conditional prior probabilities for any gene depends only on the effects of neighboring genes. The mean of the conditional distribution is the mean of the effects of neighboring genes, and its variance decrease with edge degree. A more general approach may be obtained by assigning a multivariate normal prior on gene effects with a general (co)variance structure derived from the network, in the form of a kernel.

The key element in this approach is the form of Σ, which is a function of a gene network (or graph). Figure 1 illustrates how to express an undirected unweighted graph into an adjacency matrix, E, whose entries e_ij are equal to 1 if genes (or nodes) i and j have a common edge connecting them (e.g. indicating (co)expression), or 0 otherwise. As an example, in Figure 1a, genes A and B are connected with an edge, and D is an isolated gene. Kernel candidates are functions of E, such as the exponential exp(−d_ij), or the inverse 1/d_ij, where d_ij is the length of the path between nodes i and j, with an infinite distance between no connected genes (e.g. genes A and D in Figure 1a). A particular function of E is the diffusion kernel (Kondor and Lafferty, Reference Kondor and Lafferty2002), which has been previously used in microarrays in a classification problem with a support vector machine approach (Vert and Kanehisa, Reference Vert and Kanehisa2003). The diffusion kernel matrix has larger entry values between nodes more closely connected and smaller values between two nodes with a longer path, with the strength of the correlations controlled by a bandwidth parameter h. On Figure 1, two different values of the bandwidth parameter were used to construct diffusion kernels based on the simple network depicted (h = 0.2 and 0.5 for Panels c and d, respectively).

Figure 1

Example of an undirected unweighted graph representing a 5-node network (panel a), and its corresponding adjacent matrix (panel b), and two possible diffusion kernels (panels c and d).

As an illustration of the potential of using prior biological knowledge in the analysis of microarray data, we consider a simple example with the data set and network used by Wei and Pan (Reference Wei and Pan2008). Briefly, log-ratio expression of 4609 genes of S. cerevisiae was measured in absence and presence of a TF, using ChIP-on-chip technology. A weighted average of three independent replicates in each group (i.e. with and without TF) were analyzed, incorporating or not the information provided by a probabilistic functional gene network obtained from an independent study (Lee et al., Reference Lee, Date, Adai and Marcotte2004).

The following model was used to incorporate the gene network information:

where y_k is the average (over the three replicates) of the log-ratio intensity values for gene k, m is an overall mean, g_k is the differential expression (fold change in the log scale) relative to gene k, and e_k is a residual term associated with the observation y_k. It is assumed that , and that g = [g ₁, g ₂,…,g _K]′∼N(0,K_h). A diffusion kernel with bandwidth parameter h = 0.2 (K_0.2) was then generated from the network. A model that ignores the gene network information was also considered, in which case g = [g ₁, g ₂,…,g _K]′∼N(0,Iσ_g ²), i.e. . These models were evaluated in terms of predictive ability using a leave-one-out cross-validation (LOOCV) approach. In the model with the gene network (K_0.2 model), prediction of ratio values were obtained by , whereas in the model ignoring the gene network, because of the confounding between g_k and e_k, predictions were equal to for all k.

With the K_0.2 model, the estimate of variance due to the gene structure was 0.006, whereas the residual variance was 0.045, such that the gene structure explained about 12% of the total variability. When ignoring the network, however, it was not possible to separate both variances, since the data set does not have replications. The network allowed not only the estimation of variances, but also the prediction of gene expression differences (fold changes). The correlation between predicted and observed log-ratio values was 0.76, indicating that the network produced considerable information borrowing across genes. The mean squared error of the prediction (MSEP) evaluated via LOOCV was smaller in K_0.2 compared with the model ignoring the network, being 0.0503 and 0.0545, respectively (i.e. an improvement of 7.6%).

For successful modeling of covariance structures for microarray data using prior biological knowledge, the choice of an adequate kernel is critical. Building adequate kernels requires two main components: (i) an informative input set and (ii) an appropriate function to generate the kernel from the adjacent matrix. Informative input sets can be gathered from trustable and validated biological knowledge from public databases, whereas an adequate kernel can be obtained, e.g., by comparing candidate functions using cross-validation techniques. In this case, the incorporation of prior knowledge into the statistical analysis of microarray data should yield lower MSEP, i.e. more accurate and or less biased predictions.