It has unambiguously been shown that genetic, environmental, demographic, and technical factors may have substantial effects on gene expression levels. In addition to the measured variable (s) of interest, there will tend to be sources of signal due to factors that are unknown, unmeasured, or too complicated to capture through simple models. We show that failing to incorporate these sources of heterogeneity into an analysis can have widespread and detrimental effects on the study. Not only can this reduce power or induce unwanted dependence across genes, but it can also introduce sources of spurious signal to many genes. This phenomenon is true even for well-designed, randomized studies. We introduce “surrogate variable analysis” (SVA) to overcome the problems caused by heterogeneity in expression studies. SVA can be applied in conjunction with standard analysis techniques to accurately capture the relationship between expression and any modeled variables of interest. We apply SVA to disease class, time course, and genetics of gene expression studies. We show that SVA increases the biological accuracy and reproducibility of analyses in genome-wide expression studies.
Large-scale gene expression studies allow one to characterize transcriptional variation with respect to measured variables of interest, such as differing environments, treatments, time points, phenotypes, or clinical outcomes. However, a number of unmeasured or unmodeled factors may also influence the expression of any particular gene. Besides inducing widespread dependence in measurements across genes [1,2], these influential factors create additional sources of differential expression, which, unlike gene-specific fluctuations, represent common sources of variation in gene expression that can be observed among multiple genes. We call “primary measured variables” (or primary variables) those variables that are explicitly modeled in the analysis of an expression study. These variables may or may not be associated with any given gene' s expression variation. We classify all the remaining sources of expression variation into three basic types. “Unmodeled factors” are sources of variation explained by measured variables, but are not explicitly included in the statistical model (e. g. , because their relationship to expression is intractable or the relevant measured variables were excluded because of sample size restrictions). “Unmeasured factors” are sources of expression variation that are not measured in the course of the study, so we also call these unmodeled factors. Finally, “gene-specific noise” refers to random fluctuations in gene expression independently realized from gene to gene. As a simple example meant only for illustrative purposes, consider a human expression study where disease state on a particular tissue type is the primary variable. Suppose that in addition to changes in expression being associated with disease state, the age of the individuals also has a substantial influence on expression. Thus, some genes exhibit differential expression with respect to disease state, some with respect to age, and some with respect to both. If age is not included in the model when identifying differential expression with respect to disease state, we show that this may (a) induce extra variability in the expression levels due to the effect of age, decreasing our power to detect associations with disease state, (b) introduce spurious signal due to the fact that the effect of age on expression may be confounded with disease state, or (c) induce long-range dependence in the apparent “noise” of the expression data, complicating any assessment of statistical significance for differential expression. In practice, even if age were known, it may be one of dozens of available measured factors, making it statistically intractable to determine which to include in the model. Furthermore, even measured factors such as age may act on distinct sets of genes in different ways, or may interact with an unobserved factor, making the effect of age on expression difficult to model. “Expression heterogeneity” (EH) is used here to describe patterns of variation due to any unmodeled factor. Major sources of expression variation are due to technical [3,4], environmental [5,6], demographic [7,8], or genetic [9–11] factors. It is well known that sources of variation due to experimental design or large-scale systematic sources of signal may be present in expression data [3,4, 12,13], sometimes even after normalization has been applied [14]. Genetic factors can also have a large-scale impact on gene expression levels. Specific genetic loci have been shown to influence the expression of hundreds or thousands of genes in several organisms [10,11,15]. Expression heterogeneity is particularly pronounced in human expression data, especially in the study of complex systems, such as cancer or responses to stress [16–18]. Recently, Lamb et al. proposed the “Connectivity Map” for identifying functional connections between cancer subtypes, genetic background, and drug action [19]. Lamb et al. noted EH (e. g. , due to cell type and batch effects) presented a major hurdle for extracting relevant biological signal from the Connectivity Map. In each of these studies, expression variation with respect to one or at most a handful of variables is explored. However, it is likely that in each study multiple sources of EH will act on distinct, but possibly overlapping, sets of genes. Normalization techniques are commonly used to detect and adjust for systematic expression variation due to well-characterized laboratory and technical sources [12,13,20]. However, to date there has been no approach for identifying and accounting for all sources of systematic expression variation, including variation due to unmeasured or unmodeled factors of both biological and technical sources. We show here that biological sources of variation not modeled in the analysis can be just as problematic as technical sources of variation. Here, we introduce “surrogate variable analysis” (SVA) to identify, estimate, and utilize the components of EH. Figure 1 shows the effects of failing to account for unmodeled factors in a differential expression analysis, and the potential benefits of the SVA approach. EH causes drastic increases in the variability of the ranking of genes for differential expression (Figure 1A), distorts the null distribution potentially causing highly conservative or anticonservative significance estimates (Figure 1B), and reduces the power to distinguish true associations between a measured variable of interest and gene expression (Figure 1C). However, employing SVA in these studies produces operating characteristics nearly equivalent to what one would obtain with no EH at all. We apply SVA to three distinct expression studies [7,21,22], where each study contains clear patterns of EH (Figure S1). These studies represent major classes of gene expression studies performed in practice: genetic dissection of expression variation, differential expression analysis between disease classes, and differential expression over time. We show that SVA is able to accurately identify and estimate the impact of unmodeled factors in each type of study, using only the expression data itself. We further show that SVA improves accuracy and consistency in detecting differential expression. SVA orders the significant gene lists to more accurately and reproducibly reflect the ordering of the genes with respect to their true differential expression signal. SVA is particularly useful in producing reproducible results in microarray studies, because adjusting for surrogate variables reduces differential expression due to sources other than the primary variables. These results indicate that EH is prevalent across a range of studies and that SVA can be used to capture and account for these patterns to improve the characterization of biological signal in expression analyses.
We have developed an approach called surrogate variable analysis that appropriately borrows information across genes to estimate the large-scale effects of all unmodeled factors directly from the expression data. Figure 2A shows a simulated example of EH. The primary variable distinguishes the first ten arrays from the last ten (Figure 2B); however, the unmodeled factor may have a variety of effects on expression (Figure 2C). The SVA approach flexibly captures signatures of EH, including highly irregular patterns not following any simple model, by estimating the signatures of EH in the expression data themselves rather than attempting to estimate specific unmodeled factors such as age or gender. After the surrogate variables are constructed, they are then incorporated into any subsequent analysis as covariates in the usual way. The SVA algorithm, described in mathematical detail in Materials and Methods, can conceptually be broken down into four basic steps: (Step 1) Remove the signal due to the primary variable (s) of interest to obtain a residual expression matrix. Apply a decomposition to the residual expression matrix to identify signatures of EH in terms of an orthogonal basis of singular vectors that completely reproduces these signatures. Use a statistical test to determine the singular vectors that represent significantly more variation than would be expected by chance. (Step 2) Identify the subset of genes driving each orthogonal signature of EH through a significance analysis of associations between the genes and the EH signatures on the residual expression matrix. (Step 3) For each subset of genes, build a surrogate variable based on the full EH signature of that subset in the original expression data. (Step 4) Include all significant surrogate variables as covariates in subsequent regression analyses, allowing for gene-specific coefficients for each surrogate variable. The four-step procedure is necessary both to ensure that the surrogate variables indeed estimate EH and not the signal from the primary variable (Step 1), to ensure an accurate estimate of each surrogate variable by identifying the specific subset of genes driving each EH signature (Step 2), to allow for correlation between the primary variable and the surrogate variables by building the surrogate variables on the original expression data (Step 3), and to take into account the fact that a surrogate variable may have a different effect on each gene (Step 4). The third and fourth steps are particularly important for maintaining unbiased significance with SVA, as demonstrated below. The overall goal of SVA is to provide a more accurate and reproducible parsing of signal and noise in the analysis of an expression study when EH is present. One way in which signal is commonly quantified is through a significance analysis [23]. The most basic definition of a significance analysis being performed “correctly” is if the null distribution is calculated properly [24]. A straightforward means for determining whether this is true is to assess whether the p-values corresponding to true null hypotheses are Uniformly distributed between zero and one. Indeed, p-values are specifically defined so that those corresponding to true null hypotheses have a Uniform (0,1) distribution if and only if the null distribution has been correctly calculated [25]. Throughout this paper, we examine the distribution of p-values from null genes to determine whether various procedures are able to recover the correct null distribution in the presence of EH. To assess statistically any deviations from the Uniform distribution for the null p-values, we apply a nested Kolmogorov-Smirnov test that is robust to chance fluctuations that may be present in a single simulated dataset (see Text S1). We performed a simulation study to investigate the properties of SVA with respect to large-scale significance testing. Specifically, we show that the SVA algorithm (a) accurately estimates signatures of expression heterogeneity, (b) corrects the null distribution of p-values from multiple hypothesis tests, (c) improves estimation of the false discovery rate (FDR) [23,26], and (d) is robust to confounding between the primary variable and surrogate variables. The primary variable for our simulation was a binary variable indicating two disease classes. We simulated 1,000 expression studies, drawn from the same hypothetical population. For each study, we simulated expression for 1,000 genes on 20 arrays divided between the two disease states. The first 300 genes were simulated to be differentially expressed between disease states and genes 200–500 were affected by an independent unobserved factor to simulate a randomized study (Materials and Methods). Several recent studies have carried out the genetic dissection of expression variation at the genome-wide level [10,11,15]. Brem et al. [10,21] measured expression genome wide in 112 segregants of a cross between two isogenic strains of yeast. They also obtained genotypes for each segregant at markers covering 99% of the genome (Materials and Methods). It was shown that many gene expression traits are cis-linking, i. e. , the quantitative trait locus (QTL) linkage peak coincided with the physical location of the open reading frame for the expression trait [36]. At the same time, it was also shown that a number of gene expression traits are trans-linking, with linkage peaks at loci distant from the physical location of their open reading frames. In particular, several “pivotal” loci each appear to influence the expression of hundreds or even thousands of gene expression traits. Similar highly influential loci have been observed in other organisms [11,15]. These pivotal loci act as a major source of EH, regardless of whether genotypes have been measured in an expression study. As proof of concept, the Brem et al. [10,21] dataset was used to show that well-defined EH exists in actual studies and that SVA can properly capture and incorporate this EH structure into the statistical analysis of measured variables of interest. First, we analyzed the full dataset to identify the expression traits under the influence of these pivotal trans-acting loci, as well as the patterns of EH induced by these loci. Then we applied SVA to only the expression data, ignoring the genotype data to identify relevant surrogate variables capturing EH. Linkage analysis was performed again including the surrogate variables as covariates, showing that the effects from the pivotal loci are now negligible. In other words, SVA was able to capture and remove the effects of these few pivotal loci without the need for genotypes. A number of expression traits have significant trans-linking eQTL mapping to pivotal loci on Chromosomes II, III, VIII, XII, XIV, and XV (Figure 3A). In the SVA-adjusted analysis, the majority of the trans-linkages to the pivotal loci have been eliminated (Figure 3B). The pervasive trans-linkage signal mapping to the pivotal loci can be viewed as global expression heterogeneity. The reduction in trans-linkage to these loci in the SVA-adjusted significance analysis indicates that SVA effectively captures genetic EH. Pivotal trans-linkage signals indicate large-scale effects of a few loci. However, subtle and potentially more interesting cis-linkage may be lost in the presence of substantial genetic heterogeneity. To assess the impact of SVA on power to detect cis-linkage, we calculated linkage p-values only for markers located within three centimorgans of the open reading frame of each trait. On chromosomes without a pivotal locus (Chromosomes I, IV, V, VI, VII, IX, X, XI, and XIII) the SVA-adjusted analysis finds substantially more cis-linkage signal. At an FDR cutoff of 0. 05, the adjusted analysis finds 1,894 significant cis-linkages, compared with 1,604 for the unadjusted analysis. This increase is consistent across a range of FDR cutoffs (Table 1) and illustrates the potential increase in power obtained from applying SVA. We applied the SVA approach to two human studies [7,22], representing the two common human study designs: disease state and timecourse.
Expression heterogeneity due to technical, genetic, environmental, or demographic variables is common in gene expression studies. Here we have introduced a new method, SVA, for identifying, estimating, and incorporating sources of EH in an expression analysis. SVA uses the expression data itself to identify groups of genes affected by each unobserved factor and estimates the factor based on the expression of those genes. Simulations show that SVA accurately detects expression heterogeneity and improves significance analyses. We performed SVA on experiments involving recombinant inbred lines, individuals of varying disease state, and expression measured over time to illustrate the broad range of studies on which SVA can be applied. One advantage of the SVA approach is the ability to disentangle correlated and overlapping differential expression signals. This approach may be particularly useful in clinical studies, where a large number of clinical variables may have a complicated joint impact on expression. We have implemented SVA in an open source package available for downloading at http: //www. genomine. org/sva/.
Three publicly available datasets were employed to represent a broad range of gene expression studies performed in practice. The first dataset consists of gene expression measurements for 6,216 genes in 112 segregants of a cross between two isogenic strains of yeast, as well as genotypes across 3,312 markers [10,21]. The second dataset consists of gene expression for 3,226 genes in seven BRCA1 and eight BRCA2 mutation–positive tumor samples [22]; several genes with apparent outliers were removed as described [23] for a total of 3,170 genes. The third dataset consists of gene expression measurements in kidney samples from normal kidney tissue obtained at nephrectomy from 133 patients [7]; the 34,061 genes analyzed in [8] were also analyzed here. Seventy-four of the tissue samples were obtained from the cortex and 59 from the medulla. Details of the protocol for each study appear in the corresponding references. All expression data were analyzed on the log scale. The SVA algorithm identified 14 significant surrogate variables from the expression data. We performed both an unadjusted and an SVA-adjusted linkage analysis for each expression trait. In the unadjusted analysis, linkage p-values were calculated from an F-test comparing an additive genetic model to the null model of no genetic association [42]. SVA-adjusted p-values were calculated from an F-test comparing the full model of genetic association and the null model of no association, both models including all significant surrogate variables as additive terms. For each study, we simulated expression for 1,000 genes on 20 arrays divided between the two disease states. For simplicity, the expression measurements for each gene were drawn from a normal distribution with mean zero and variance one. We simulated expression heterogeneity with a dichotomous unmodeled factor independent of the disease state. The mean differences between disease states and states of the unmodeled factor were drawn from two independent normal distributions. For the real data example, we calculated the residuals from the regression of BRCA tumor type on expression for the Hedenfalk data [22]. Then, for each simulated study, we independently permuted each row of the expression data to create a matrix of residuals. To this matrix, we added signal, as in the case of the purely simulated data. The simulation studies were based on data generated using the R programming language [43]. All differential expression analyses were performed by a t-test based on standard linear regression. The genes were ranked for relative significance by the absolute values of their t-statistics. Differential expression was calculated using a t-test based on standard linear regression for the disease class data. The method of Storey et al. [8] was applied for the time-course data. q-Values were estimated using previously described methodology [23]. Let Xmxn = (x1, . . , xm) T be the normalized m × n expression matrix with n arrays for m genes, where xi = (xi1, . . , xin) T is the vector of normalized expression for gene i. Let y = (y1, . . , yn) T be a vector of length n representing the primary variable of interest. Without loss of generality model xij = μi + fi (yj) + eij, where μi is the baseline level of expression, fi (yj) = E (xij | yj) − μi gives the relationship between measured variable of interest and gene i, and eij is random noise with mean zero. As a simple example, for a dichotomous outcomes yj ∈ {−1,1} we would employ the linear model xij = μi + βi yj +eij and estimate μi and βi by least squares. We could then perform a standard test of whether βi = 0 or not for each gene. This hypothesis test is exactly equivalent to performing a test of differential expression between the two classes. Suppose in a microarray study there are L biologically meaningful unmodeled factors, such as age, environmental exposure, genotype, etc. Let gℓ = (gℓ1, . . . , gℓn) be an arbitrarily complicated function of the ℓth factor across all n arrays, for ℓ=1,2, . . . , L. Therefore, we can now model the expression for gene i on array j as xij = μi + fi (yj) +, where γℓi is a gene-specific coefficient for the ℓth unmodeled factor. If unmodeled factor ℓ does not influence the expression of gene i, then γℓi = 0. The fact that we employ an additive model is actually quite general: it has been shown that even complicated nonlinear functions of factors can be represented in an additive fashion for a reasonable choice of a nonlinear basis [44]; we simply define the gℓ to be as nonlinear as necessary and make L as large as necessary to best fit the additive effect. Since there are n arrays, each gene' s expression can be modeled by at most n linearly independent factors, and hence any dependence structure between genes can be represented using L ≤ n vectors in this additive fashion. Due to this formulation, the inter-gene dependent eij have now been replaced with, where is the true gene-specific noise, now sufficiently independent across genes. In other words, we have broken the error eij into two terms, one that represents dependent variation across genes due to unmodeled factors, , and one that represents gene-specific independent fluctuations in expression. It is not possible in general to directly estimate the unmodeled gℓ, and SVA does not attempt to do so. The key observation is to note that for L vectors in n space, it is possible to identify an orthogonal set of vectors hk, k = 1, . . . , K (K≤L) that spans the same linear space as the gℓ In other words, for any set of vectors gℓ and coefficients γℓi, it is possible to identify mutually orthogonal vectors hk and coefficients λki such that and Therefore, we do not need to estimate the specific variables gℓ. We only need to estimate the linear combination, so we can choose a set of vectors that spans the same space but is statistically tractable. Here we choose the set of K orthogonal vectors (denoted by the hk) to be those that are the right non-zero singular vectors provided by the singular value decomposition of the m × n matrix with (i, j) entry. This justifies the use of the singular value decomposition to identify orthogonal signatures of expression heterogeneity for surrogate variable estimates. We call these h1, h2, . . . , hK the “surrogate variables. ” An intuitive question that arises from an inspection of this formulation is about the model assumptions of the gℓj. Whereas the term fi (yj) is a model of the measured variable, yj, it is not generally possible to analogously formulate gℓj as a function of a well-defined, measured variable. Since we estimate the outcomes directly from the expression data (as), it is not necessary to determine a model of the gℓj in terms of a biologically meaningful variable. Thus, we can bypass the need to know what the most relevant model of a measured variable is for gℓj for the purposes of estimating the EH. The goal of the SVA algorithm is therefore to identify and estimate the surrogate variables, hk, = (hk1, . . . , hkn) T, based on certain consistent patterns of expression variation. Methods for empirically identifying [37] and estimating [40] expression trends or clusters have previously been developed. However, care must be taken when estimating expression trends for use in subsequent analyses of measured variables of interest. Specifically, surrogate variables must represent signal due to sources other than the primary variable and allow for potential overlap with the primary variable. The SVA algorithm is designed to estimate surrogate variables that meet both requirements. We assume that n < m and, for simplicity, that there is only a single primary variable; the extension to multiple primary variables simply requires one to include all of them in the model fit occurring in each Step 1 below. The algorithm is decomposed into two parts: detection of unmodeled factors and construction of surrogate variables. The basic form of the first algorithm has been proposed previously [27]. The second algorithm has been proposed and justified in this manuscript Algorithm to detect unmodeled factors. 1. Form estimates and by fitting the model xij = μi + fi (yj) + eij, and calculate the residuals rij = xij − (yj) to remove the effect of the primary variable on expression. Form the m × n residual matrix R, where the (i, j) element of R is rij. 2. Calculate the singular value decomposition of the residual expression matrix R = UDVT. 3. Let dℓ be the ℓth eigenvalue, which is the ℓth diagonal element of D, for ℓ=1, . . . , n. If df is the degrees of freedom of the model fit μ̂i + f̂i (yj), then by construction the last df eigenvalues are exactly zero and we remove them from consideration. For eigengene k=1, . . . , n-df set the observed statistic to be which is the variance explained by the kth eigengene. 4. Form a matrix R* by permuting each row of R independently to remove any structure in the matrix. Denote the (i, j) entry of R* by. 5. Fit the model and calculate the residuals to form the m × n model-subtracted null matrix R0. 6. Calculate the singular value decomposition of the centered and permuted expression matrix R0 = U0D0. 7. For eigengene k form a null statistic as above, where d0ℓ is the ℓth diagonal element of D0. 8. Repeat steps 4−7 a total of B times to obtain null statistics for b = 1, . . . , B and k = 1, . . . , n-df. 9. Compute the p-value for eigengene k as: Since eigengene k should be significant whenever eigengene k′ is (where k′>k), we conservatively force monotonicity among the p-values. Thus, set pk = max (pk−1, pk) for k = 2, . . . , n-df. 10. For a user-chosen significance level 0≤α≤1, call eigengene k a significant signature of residual EH if pk ≤ α. Algorithm to construct surrogate variables. 1. Form estimates and by fitting the model xij. = μi + fi (yj) + eij, and calculate the residuals rij = xij − (yj) to remove the effect of the primary variable on expression. Form the m × n residual matrix R, where the (i, j) element of R is rij. 2. Calculate the singular value decomposition of the residual expression matrix R = UDVT. Let ek = (ek1, . . . , ekn) T be the kth column of V (for k=1, . . . , n). These ek are the residual eigengenes and represent orthogonal residual EH signals independent of the signal due to the primary variable. 3. Set to the number of significant eigengenes found by the above algorithm. Note that “significant” means that the eigengene represents a greater proportion of variation than expected by chance. For each significant eigengene ek k=1, . . . , . 4. Regress ek on the xi (i = 1, . . . , m) and calculate a p-value testing for an association between the residual eigengene and each gene' s expression. This p-value measures the strength of association between the residual eigengene ek and the expression for gene i. 5. Let π0 be the proportion of genes with expression not truly associated with ek; form an estimate, as described previously [23] and estimate the number of genes associated with the residual eigengene by. Let be the indices of the genes with smallest p-values from this test. 6. Form the × n reduced expression matrix. Since is an estimate of the number of genes associated with residual eigengene k, the reduced expression matrix represents the expression of those genes estimated to contain the EH signature represented by some hk as described above. As was done for R, calculate the eigengenes of Xr, and denote these by for j=1, . . . , n. 7. Let j* = argmax1≤j≤n cor and set. In other words, set the estimate of the surrogate variable to be the eigengene of the reduced matrix most correlated with the corresponding residual eigengene. Since the reduced matrix is enriched for genes associated with this residual eigengene, this is a principled choice for the estimated surrogate variable that allows for correlation with the primary variable. 8. In any subsequent analysis, employ the model xij = μi + fi (yj) +, which serves as an estimate of the ideal model xij = μi + fi (yj) +. The singular value decomposition is employed in these SVA algorithms. It may be possible to utilize other decomposition methods, but since the singular value decomposition provides uncorrelated variables that decompose the data in an additive linear fashion with the goal of minimizing the sum of squares, we found this to be the most appropriate decomposition. If the primary variables are modeled for data that are not continuous, then it may make sense to decompose the variation with respect to whatever model-fitting criteria will be employed SVA has been made freely available as an R package at http: //www. genomine. org/sva/.
In scientific and medical studies, great care must be taken when collecting data to understand the relationship between two variables, such as a drug and its effect on a disease. In any given study there will be many other variables at play, such as the effects of age and sex on the disease. We show that in studies where the expression levels of thousands of genes are measured at once, these issues become surprisingly critical. Due to the complexity of our genomes, environment, and demographic features, there are many sources of variation when analyzing gene expression levels. In any given study, it is impossible to measure every single variable that may be influencing how our genes are expressed. Despite this, we show that by considering all expression levels simultaneously, one can actually recover the effects of these important missed variables and essentially produce an analysis as if all relevant variables were included. As opposed to traditional studies, the massive amount of data available in this setting is what makes the method, called surrogate variable analysis, possible. We hypothesize that surrogate variable analysis will be useful in many large-scale gene expression studies.
saccharomyces
homo (human)
computational biology
Capturing Heterogeneity in Gene Expression Studies by Surrogate Variable Analysis