A factor model to analyze heterogeneity in gene expression

Background Microarray technology allows the simultaneous analysis of thousands of genes within a single experiment. Significance analyses of transcriptomic data ignore the gene dependence structure. This leads to correlation among test statistics which affects a strong control of the false discovery proportion. A recent method called FAMT allows capturing the gene dependence into factors in order to improve high-dimensional multiple testing procedures. In the subsequent analyses aiming at a functional characterization of the differentially expressed genes, our study shows how these factors can be used both to identify the components of expression heterogeneity and to give more insight into the underlying biological processes. Results The use of factors to characterize simple patterns of heterogeneity is first demonstrated on illustrative gene expression data sets. An expression data set primarily generated to map QTL for fatness in chickens is then analyzed. Contrarily to the analysis based on the raw data, a relevant functional information about a QTL region is revealed by factor-adjustment of the gene expressions. Additionally, the interpretation of the independent factors regarding known information about both experimental design and genes shows that some factors may have different and complex origins. Conclusions As biological information and technological biases are identified in what was before simply considered as statistical noise, analyzing heterogeneity in gene expression yields a new point of view on transcriptomic data.


Background
Microarray technology allows the analysis of expression levels for thousands of genes simultaneously and is a powerful tool to characterize mRNA level variation due to measured variables of interest (various phenotypes, treatments...). Typical approaches to find significant relationships between gene expressions and experimental conditions ignore the correlations among expression profiles and functional categories [1]. This dependence structure leads to correlation among test statistics which affects a strong control of the actual proportion of false discoveries [2]. Indeed, a number of unmeasured or unmodeled factors independent of the variables of interest may influence the expression of any particular gene [3,4]. These factors may induce extra variability in the expression levels and decrease the power to detect links with the variables of interest.
Recently, several works have introduced models for the common information shared by all the genes. Especially Friguet et al [4] propose to model this sharing of information by a factor analysis structure in a method called Factor Analysis for Multiple Testing (FAMT). The estimated factors in the model capture components of the expression heterogeneity. As well, Storey et al [3] introduce Surrogate Variable Analysis (SVA) to identify and estimate these extra sources of variation. The factors in FAMT and the surrogate variables in SVA are similarly designed to model dependence among tests by a linear kernel but they are estimated differently. Contrarily to the SVA model, independence between the factors and the experimental conditions of interest is explicitly assumed in FAMT in order to separate clearly the effects of the experimental conditions on the gene expressions and the nuisance variability due to unmodeled technological effects and other known or unknown effects that could be uncontrolled in the experimental design.
The major sources of expression variation are then assumed to be the experimental conditions of interest, but also gene dependence and uncontrolled factors in the experimental design. Indeed, even after normalization, variation due to the experimental design still exists in expression data. The factors extracted in the residual part of the regression models explaining the gene expressions by the experimental conditions of interest are therefore analyzed to give more insight both on expression heterogeneity among sampling units and the contribution of some biological processes to gene dependence. First, factors are extracted from illustrative expression data sets with simple patterns of expression heterogeneity in order to show how they can straightforward be related to sources of heterogeneity. Henceforth, the same factor model approach is used to analyze an expression data set initially generated to map quantitative trait loci (QTL) for abdominal fatness (AF) in chickens, especially on chromosome 5 (GGA5) [5]. This data set concerns hepatic transcriptome profiles for 11213 genes of 45 half sib male chickens generated from a same sire. This sire was generated by successive inter-crossing of two experimental chicken lines divergently selected on AF and was known to be heterozygous for an AF QTL on the GGA5 chromosome around 175 cM (For more details, see [5]). The 45 half sib chickens show therefore variation on AF. According to the polygenic effect model of quantitative traits, this variation is probably due to multiple mutations and biological processes.
Two lists of genes significantly correlated to the AF trait are first generated using the raw and the factor-adjusted expression dataset. Then, the relevance of the two gene lists to characterize functionally fatness variation in the family are compared, regarding the frequencies of biological processes related to the AF trait in their functional annotations. Factor-adjusted expression data is finally used to identify a gene whose expression is controlled by the AF QTL region.
Furthermore, the extracted factors are interpreted using external information on the experimental design such as the hatch, dam and body weight and also gene information such as functional categories, oligonucleotide size and location on the microarray. It is deduced that some factors may have different and complex origins, which confirms the importance of taking into account these extra sources of variability to be more relevant in the transcriptomic analyses.

Illustrative Examples
Similarly to Storey et al (2007) [3], three simple situations of heterogeneity are considered. For each one, independent expressions for 1000 genes on 20 arrays are simulated according to a standard normal distribution. The sample is split into two equal groups and a constant is added on the first 100 gene expressions to mimic a differential expression between these two groups.

Case 1: One independent variable affecting all genes
All genes are affected by an independent grouping variable marked by colors red and green on Figure 1A. A single factor is extracted by FAMT. Figure 2A helps interpreting this factor and shows that it clearly discriminates the two colored groups of individuals (P-value ≤ 2.2 × 10 -16 ). This shows a high association between the factor and the independent grouping variable. The genes representation does not show any particular structure. In this simple case the factor estimated by FAMT can therefore be easily interpreted through the individuals representation.

Case 2: One independent variable affecting a set of genes
Only genes 70-170 are affected by an independent grouping variable marked by colors red and green on Figure 1B. A single factor is also found using FAMT. As shown on Figure 2B, the factor discriminates the two groups of individuals (P-value ≤ 2.2 × 10 -16 ) and the two groups of genes (P-value ≤ 2.2 × 10 -16 ). In this case, the estimated factor can be interpreted through the individuals and genes representations.

Case 3: Two independent variables affecting two different sets of genes
Gene sets 70-170 and 171-271 are each affected by an independent grouping variable marked respectively by colors red and green and by colors orange and blue as illustrated by Figure 1C. Two factors are identified by FAMT which are now interpreted regarding the two external sources of heterogeneity ( Figure 2C). The redgreen variable seems to be highly associated with the first axis (P-value ≤ 2.2 × 10 -16 in both representation). On the contrary, the orange-blue variable is not associated with this axis considering a significance level of 0.05 (P-value = 0.7933 for the individuals representation, p-value = 0.1109 for the genes representation). The same strategy is implemented for the second factor. The red-green variable appears to be not associated with this factor (P-value = 0.7949 for the individuals representation, p-value = 0.1926 for the genes representation) whereas the orangeblue variable is highly associated (P-value ≤ 2.2 × 10 -16 in both representations). In this case, each of the two estimated factors can be explained by one of the two independent grouping variables.

Analysis of the AF expression data set Classical approach
Examination of the Pearson coefficient correlation between hepatic transcript levels and AF trait shows that 287 genes are significantly correlated considering a significance threshold of 0.05 without any correction for multiple tests. This low amount of differentially expressed genes might be explained by a poor genetic variability between individuals which are half sib offsprings and could also be due to dependence between genes that can lead to under representation of the smallest p-values [6].

Heterogeneity analysis
Minimizing the variance inflation criterion proposed by Friguet et al. [4], six factors containing a common information shared by all genes and independent from the AF trait are extracted. Subtracting the linear dependence kernel defined by these factors from the raw expression data yields the factor-adjusted expression data. The significance analysis based on these expressions results in a list of 688 gene expressions significantly correlated to the AF trait. 93% of the 287 genes found with the classical approach are included in this list. This larger number of differentially expressed genes suggests that correlation between many gene expressions and the variable of interest is under estimated due to gene dependence. Considering the Gene Ontology (GO) terms and KEGG pathways, one enriched term related to the lipid metabolism is found in the gene list resulting from factor-adjustment (688 genes) whereas none is observed in the gene list obtained using the raw expressions (287 genes). This term concerns "Steroid biosynthesis process" with 3 genes associated (Table 1). More precisely, these genes are involved in the cholesterol metabolism or in conversion of cholesterol in steroids. Several works show relationships between cholesterol metabolism and obesity [7][8][9]. This result shows that the genes found after factoradjustment are more related to the fatness trait. Furthermore, the impact of factor-adjustment is shown in Figure  3, where a principal component analysis (PCA) generated with the 688 factor-adjusted transcript levels of correlated genes ( Figure 3B) separates much more fat and lean chickens than the same PCA generated with the raw expressions of the same 688 genes ( Figure 3A). This observation displays that factor-adjustment has cleaned up the data from dependence, which highlights masked relationships with the AF trait.
We focus on one of the 3 genes involved in the "Steroid biosynthesis process", DHCR7, which is only observed in the list of 688 genes and known for encoding the last enzyme involved in the cholesterol synthesis. As shown in Figure 4, the analysis of the factor-adjusted expressions for this gene highlights an eQTL (P-value < 0.05) colocalizing with the AF QTL previously observed [5]. The same LRT curve based on the raw expressions does not point out any eQTL. This result shows that the expression of this gene is controlled by a mutation in the same GGA5 AF QTL region. Because of the function of this gene related to lipid metabolism, this result suggests that this mutation could be the same as the QTL mutation for fatness phenotype. Further investigations are necessary to refine these QTL and eQTL locations.

Factor interpretation
In the present study, some external information about the experimental design and the genes is available. As we did in the simulated examples, we interpret the factors extracted from the AF expression dataset using this known information.

Using information on experimental design
The hatch, the dam and the body weight were previously measured for each bird and should be independent of the AF variation. For the body weight, the founder chicken lines were selected on AF criteria maintaining a constant body weight. The variables "hatch" and "dam" are both categorical with respectively, four and eight levels and the body weight is a continuous variable. We first focus on the "hatch" and for each factor we represent the individuals colored according to their hatch ( Figure 5). Factor 1 seems to discriminate hatches 1 and 4, factor 3 hatch 2 from the others and factor 4 hatches 2 and 3 from hatch 1. The effect of the hatch on each factor is tested and the results given in Table 2 confirm our previous observations: factor 1, 3 and 4 can be partly explained by a hatch effect (the significant test for each hatch level is given in Additional file 1). We then calculate the association for the "dam" and "body weight" with each of the six factors. Table 2 shows no effect of the dam and a high correlation between the weight and factor 2. Contrarily to the illustrative cases where each factor could be interpreted by a unique variable, the factors found here seem to have more complex origins. Indeed, three of the six factors can be interpreted by an hatch effect and another one by a body weight effect. The same analysis is now performed after adjustment of the raw expression data for hatch and body weight. Interestingly, only five factors independent of the AF trait are extracted and still a hatch effect exists but only on the first factor and a weight effect on the second factor ( Table 2). This persistence of both effects suggests that there exists an interaction involving hatch and body weight with other unmeasured and/or unknown variables. Therefore, taking into account the hatch and body weight in the statistical model seems to be not sufficient to remove a consequent part of the heterogeneity in gene expression.

Using gene information
To interpret the estimated factors in terms of gene expressions, we use known information about genes as oligonucleotide size and location on the chip: block, row and column (genes representation on each factor is given in Additional file 2). For these variables, we test their association with each factors extracted from the raw data and the hatch and body weight adjusted data. As shown in Table 2, there is a strong oligonucleotide size effect and block effect captured by almost all the factors. We exhibit also a row and column effect associated with some factors. Moreover, the genes that most contribute to the first two factors are identified (score larger than 0.8). We obtain a set of 313 genes for factor 1 and a set of 175 genes for factor 2 and which were as expected not included for 95% of them in the list of 688 genes. We perform a term enrichment test for this two sets ( Table 3). As we expect, there are essentially biological terms independent of the lipid metabolism. Factor 1 is mainly characterized by genes involved in cell division metabolism and interestingly also to pigmentation. Factor 2 is more characterized by genes involved in the nucleotide metabolism. The enriched terms found are thus not implicated in the metabolim changes induced by the AF variability. We previously highlighted a hatch and body weight effects on the factors. As in PCA, individuals and genes representa-   tions can be interpreted commonly. Hatch effect could therefore be related to the particular metabolisms characterizing factor 1.

Discussion and Conclusion
The model used in the present study assumes that the gene expressions are uncorrelated given a set of hidden variables called factors. In comparison to classical methods which do not take into account the dependence between genes, this approach provides a list of genes more correlated to the variable of interest. Moreover, factor-adjustment of the expression dataset turns out to give more insight to subsequent analyses such as QTL characterization. As a result, a gene is identified as correlated to the AF trait and related to the cholesterol metabolism having a trans-eQTL colocalizing with GGA5 AF QTL.
Because several works show a link between cholesterol and obesity, this gene could be considered as a signature of the mutation underlying this AF QTL rather than a mutation close to it. This result provides functional hypothesis about genes whose expression could be impacted by the QTL of interest. Factor analysis was introduced in the psychometric field in 1904 by Spearman [10] in order to extract the common factors in intelligence and personality. In this particular domain, the individuals are explained by their responses to different subsets of tests. The method usually furnished at least five factors which were interpreted as follows: neuroticism, extraversion, conscientiousness, agreeableness and openness to ideas. In our study, the factors were found using an EM algorithm presented by [4]. Our purpose was first to interpret the estimated factors and consequently to investigate which kind of information present in this factor structure could generate heterogeneity of the gene expressions. External information concerning the experimental design and functional annotations of the genes were used to analyse the factors. It is deduced that some factors seem to have a complex explanation with at least 2 variables associated to them. For factor 1, the individuals variability independent of the trait of interest is for instance shown to be related to the hatch. Enrichment tests also give a characterization of this factor by specific metabolisms.
To remove expression heterogeneity from the data for the subsequent statistical analyses, the basic idea consists in adjusting the raw expression data from the common factor structure. As we extract uncontrolled effects and technological biases from what was before simply consid- ered as statistical noise, analyzing heterogeneity in gene expression yields a new point of view on transcriptomic data. We show in this study the importance of taking into account these extra sources of variation to be more relevant in the transcriptomic analyses.

AF expression data set
The data set concerns hepatic transciptome profiles for 11213 genes of 45 half sib male chickens variable for abdominal fatness (AF). The data set was generated to map quantitative trait loci (QTL) for abdominal fatness in chickens and used in a previous study [5]. The sire of this family, generated by successive inter-crossing of two experimental chicken lines divergently selected on abdominal fatness, was known to be heterozygous for an AF QTL on the GGA5 chromosome around 175 cM. Animals, marker genotyping and transcriptome data acquisition and normalization are described in Le Mignon et al (2009) [5].

Illustrative examples
For each case, we simulated expression for 1000 genes on 20 arrays divided in two groups using the R programming language. Initially, the expression measurements for each gene were independently drawn from a standard normal distribution. The expression heterogeneity due to simple independent grouping variables was included in the simulated data set by adding a constant value for 7 random individuals for all genes (case 1) or a set of genes (case 2 and 3).

Classical expression analysis
As the variable of interest in the biological study is continuous, we calculated the Pearson correlation coefficient for each gene expression and deduced the number of genes correlated to the trait by considering the P-values under the cutoff 0.05.

Factor-analytic method Steps and algorithm
The method takes into account the impact of dependence on the multiple testing procedures for high-throughput The p-value is given for each association test (see Methods section). Considering a threshold of 1%, the significant p-values are in bold.
data. The common information shared by all the variables (i.e. gene expressions) is modeled by a factor analysis structure. Let Y (k) = (Y (1) , Y (2) ,..., Y (m) )' be a random m-vector and x (k) = (x (1) ,..., x (p) )' some explanatory variables. The conditional covariance matrix of the responses, given the explanatory variables, is represented by a factor analysis model: Σ = Ψ + BB', where Ψ is a diagonal m × m of uniquenesses and B is a m × q matrix of factor loadings.
In the above decomposition, the diagonal elements in Ψ are referred to as the specific variances of the responses and therefore BB' appears as the shared variance in the common factor structure. This factor analysis representation of the covariance is equivalent to the following mixed effects regression modeling of the data: for k = 1, where b k is the kth row of B, Z = (Z (1) , ..., Z (q) ) are latent factors supposed to concentrate the common information in the m-responses and e = (e (1) ,..., e (m) )' is a normally distributed m-vector independent of Z, with mean 0 and variance-covariance Ψ.
An EM algorithm [11] is used to estimate Ψ, B and Z. The number of factors is chosen so that the variance of the number of false discoveries is minimized. A VARI-MAX rotation is finally applied on the factors after EM estimation in order to privilege highly dispersed loadings rather than a homogeneous distribution of the loadings. Once the factor model is estimated, factor-adjusted test statistics are obtained by correction of the classical tests from the effect of the common factors. [4] show that the Enrichment tests were performed on the genes contributing the more to the construction of factor 1 and 2 using the GO BP terms and Kegg pathways. For each enriched term, the identifier (ID), the bioligical term, the size of the whole list of genes related to the term (size), the number of genes in the sub-list related to the term (count), the pvalue of the test.
resulting tests statistics are asymptotically uncorrelated, which improves the overall power of the multiple testing procedure. The algorithm is implemented in the "FAMT" R package available from CRAN. For the subsequent analyses, the raw expression data set is adjusted for the estimated independent factors, which results in the socalled factor-adjusted expression data :

Individual and variable representation
As in PCA, the data set is transformed into a new coordinate system by an orthogonal linear transformation [12]. We can represent the individuals and variables graph through B, the matrix of factor loadings and Z, the matrix of estimated factors. Those two representations are related by a transition formula [12], which enables their simultaneous interpretation. Moreover, each factor can be related to external information which may be available in the experimental design (significance of the relationship is assessed by an analysis of variance test).

QTL and eQTL mapping
QTLMAP software based on an interval mapping method described by Elsen et al [13], was used to detect QTL affecting the AF trait and the eQTL affecting the expression of DHRC7. The statistical variable for testing the presence of one QTL (or eQTL) versus no QTL (or no eQTL) at one location was an approximate likelihood ratio test (LRT) [14]. Significance thresholds were empirically determined for AF QTL and DHCR7 eQTL from 2000 simulations. For more details, see Le Mignon et al (2009) [5].

Gene set enrichment
The enrichment of biological terms among a list of genes was assessed by the probability that an equally high or higher enrichment could be obtained by chance given the frequency of the biological terms among all the genes considered. We first implemented an R program which calculated the P value using a Fisher exact test for overrepresentation and return the enriched terms. Let denote the subset of genes related to a given metabolism in a gene set of interest. The Fisher exact test corresponds to the hypergeometric sum as follows: where . B the number of genes contained in the whole population, m the number of genes in the gene set of interest and B 0 the number of genes related to the metabolism. The functional annotations used for this program were generated as indicated in [15] are available on the website: http:// www.sigenae.org. They were obtained by a bioinformatics procedure using the Ensembl annotation source [16]. The analysis were done using the Gene ontology (GO) biological processes (BP) terms [17] and the KEGG pathways [18] with a significant threshold of 0.05. Additional material