Gene-expression measurement: variance-modeling considerations for robust data analysis

The error structure of microarray data

Measurement of variance is an indication of the reproducibility of several measurement obtained with the same experimental conditions. Many potential sources of variability exist, including noise introduced by the experimental equipment and fluorescence quality, which we categorize as gene expression–independent variability, as well as hybridization strength and RNA-input amount, which we categorize as gene expression–dependent variability. Biological variability (the changes in gene expression that are variable between controlled measurements) may contribute to both classes of variance. Both gene expression–independent and gene expression–dependent variability lie beyond the control of the experimenter, and many of these variables can influence individual genes differently. It has been widely observed that the degree of measurement uncertainty is dependent on the amount of gene expression and, in many cases, the variance decreases as a fraction of the total signal as expression increases. At greater signal intensities, both classes of variability are small relative to the signal. This error structure of high-throughput data can itself be a reliable measure of noise and may be used to devise tests with better accuracy. For determining whether a particular gene measured on an array has different expression in two physiological conditions, it is necessary to account for uncontrollable sources of variation.

To offset the typically poor estimation of variance and to increase statistical power, investigators have developed methods that borrow information from genes across the array. Distribution parameters or pseudo-replicates can be shared for genes, with the assumption that different genes may have similar patterns of variability. This strategy is largely agreed on by many groups and has been implemented in many software packages, including but not limited to SAM (significance analysis of microarrays)⁵, Cyber-T⁶, Limma (linear models for microarray data)⁵ and VAMPIRE (variance-modeled posterior inference with regional exponentials)^8,9. Detailed descriptions of the statistical modeling that underlies these analysis programs have been provided elsewhere^5–9. SAM augments estimates of the standard deviation with an estimated positive constant, so that variances closer to zero are proportionately penalized. Cyber-T uses a posterior variance estimate based on a Bayesian probability conditional on observed data, which combines sample variance and a background pooling of variances from genes with similar expression intensities contained in a window size defined by the user. Limma instead uses a posterior variance that represents a weighted combination of sample variance and a more stable grand variance estimate that is approximated from all data on the array. VAMPIRE is distinct in that it uses variance estimates derived from a global model of the intrinsic variance structure of microarray data. The global model that describes the relationship between expression intensity and variance (σ²) can be written as follows:

where A is the expression-dependent variance, B is the expression-independent variance and mu is the true expression. The Markov chain Monte Carlo (MCMC) algorithm is used for the computation of estimates for the hyper-parameters A and B. As with least-squares regression, μ is not known before the observation of array measurements; thus, the observed sample mean is often used as a surrogate. This choice, although necessary, can be problematic at low intensities at which poor signal-to-noise ratios can cause underestimation of true variability. Therefore, estimation techniques that identify parameters that are stable relative to an expression cutoff are often used. From this global variance model, estimates of variance in gene expression are extrapolated and integrated into a statistical test.

In a two-channel microarray, measurements are obtained from two RNA samples, each labeled with a different dye. These samples are subsequently hybridized to microarray probes on a single array. Typically, one RNA sample represents a control and the other represents some kind of experimental treatment. In this system, we refer to each array as a single replicate (n = 1). Parallel RNA samples from the control and experimental group can be subsequently labeled in reverse and hybridized to a second array to create a ‘dye-swap’ replicate pair. These may either be aliquots of the previously used RNA samples (for array replicates) or samples from parallel experiments (for biological replicates). The simultaneous quantification of gene expression by two RNA samples leads to covariation of the measured signals. This covariation can be modeled by the introduction of correlation coefficients for expression-dependent and expression-independent variance. A variance model of this type contains a total of six parameters. These parameters are related to the treatment condition variances (σ₁ and σ₂) and correlation (ρ) through the following equations:

where σ₁ and σ₂ are variance in the control and experimental sample, respectively, for a single probe set; A₁ and A₂ are the expression-dependent variances; B₁ and B₂ are the expression-independent variances; ρ_A is the expression-dependent correlation; ρ_B is the expression-independent correlation; and μ₁ and μ₂ are the true expression in each treatment condition. Parameters for the variance model can be computed by MCMC simulation to determine maximum-likelihood estimates for each of the six variance parameters, as achieved with the VAMPIRE microarray analysis web application. Standard statistical testing and significance thresholds that correct for multiple-parallel tests are deployed in variance-modeling strategies^6–9.

Sample size and variance

Small sample size remains a characteristic feature of most high-throughput biology experiments. This leads to unstable estimates of measurement variances after data analysis and thus leads to inaccurate detection of measured concentration differences of biomolecules. This is evident in the increasingly divergent results obtained with different statistical methods as sample size decreases. For the Student’s t-test and Welch’s t-test, small sample size is especially debilitating, as it leads to a low-powered test statistic¹⁰.

To illustrate unstable variance estimation at small sample size, we present a two-replicate subset from the Affymetrix Latin Square data set¹¹ processed along with a complete data set composed of 42 replicates (Fig. 1a). The sample variance for two replicates has wide vertical scatter for features with low expression. Other two-replicate combinations show similar poor fit. As expected, the sample variance estimates are lower with 42 replicates. The superimposed VAMPIRE plots show that these two-replicate variance estimates are more stable and accurate than standard sample variance for small sample size. In addition to modeling the expression-to-variance relationship, VAMPIRE’s model penalizes features with low expression with an additive variance component.

The importance of variance modeling of microarray data sets with fewer replicate measurements is further illustrated by biological analysis of RAW264.7 mouse macrophages treated with an inflammatory stimulus (Fig. 1b). Analysis of data from three pairs of dye-swapped samples across six time points⁹ shows that VAMPIRE identifies a nearly equal number of genes with statistically significant regulation with as little as a single pair of dye-swapped measurements as it would if six replicates were used. Cyber-T requires more replicates for identifying the same number of genes with a significant result (Fig. 1b). Methods that involve mean-based approaches, and even SAM, miss hundreds of highly significant gene-expression changes.

Sensitivity in variance modeling

The most incisive view of the importance of variance in microarray measurements has been provided by the analysis of highly precise microarrays in which a known abundance of each RNA is hybridized to the microarray (‘spike-in’ samples). A Platinum Spike data set from an Affymetrix array allows investigators to characterize in detail the repeatability of their analysis pipeline¹². In one such analysis, over 2,000 spike-in RNAs, representing expression changes from 1.2- to fourfold in two experimental conditions (A and B), are mixed with a background of over 3,500 equally abundant RNAs. To simulate an experiment with a smaller sample size, we present a three-replicate subset of the Platinum Spike as an MA plot (dye-pair intensity ratio (M) plotted by average intensity (A); Fig. 2. Here, the spikes are increasingly difficult to detect as average expression decreases. In the region of low expression, not only are the spikes embedded with background and equally spiked features, but also the upregulated and down-regulated features lie on both sides of the axis, thereby leading to potentially false-positive results (Fig. 2a). With the use of VAMPIRE at a multiple-comparison corrections Bonferroni cutoff of 5% (α_Bonf = 0.05), 77% of the spikes are detected correctly, with a false-positive rate of 2%. Cyber-T detects 71% of the spikes at the same cutoff, and the t-test detects 4% of these (data not shown).

We also present the results of an RNA spiked in equally in conditions A and B and therefore with no differences in expression (Fig. 2b). This equally spiked feature is found to be significant by the t-test but not by VAMPIRE, as VAMPIRE’s variance model systematically modulates low-expression measurements. Thus, this mRNA expression is a false-positive result incorrectly identified as being significant by t-test analysis of the microarray data set. For comparison, we present the results of a different RNA spiked in differently in conditions A and B and therefore with truly different expression (Fig. 2c). This true-positive result is considered to represent different expression by the t-test, which is correct. VAMPIRE, which results in less modulation of measurements of higher expression (and therefore measurements with more confidence), also considers this probe set to represent different expression. Given these results of spike-in experiments with well-defined constant background cRNA, variance-modeling approaches such as VAMPIRE and Cyber-T perform well in unpaired and paired experimental design settings.

Variance modeling and fold change

It is widely recognized for gene-expression data that the signal-to-noise ratio is far worse for measurements of lower intensity than for those of higher intensity. Because of this, microarray users often define an arbitrary threshold, statistical or otherwise, to remove data of low reliability before analysis. Then, an observed change of twofold or greater is deemed significant.

Although it is not a statistically robust approach, the fold-change method can give fairly reliable results, detecting many genes known to be regulated by their corresponding stimuli. Still, this approach is limited because of the need to specify arbitrary thresholds that are presumed to indicate significance. Furthermore, these thresholds are rarely, if ever, determined a priori or correlated with statistical concepts of significance. These thresholds are therefore subject to experimenter bias in identifying genes with differences in expression. To illustrate the consequences of using a fold-change approach, we present an analysis of genes that are regulated differently in adipose tissue of human patients who are insulin resistant and of a control population who are insulin sensitive^13,14 (Fig. 3).

The region that falls short of the twofold range contains a rich repertoire of genes related to the immune system established as being important in insulin resistance in human adipose tissue (Fig. 3a). Conversely, genes identified as having a difference in expression of over twofold relate to neuroactive ligand-receptor interactions or regulation of actin cytoskeleton and are not biologically important (Fig. 3a). We also present results for two pathways related to the immune system showing statistically significant enrichment: antigen processing and presentation and leukocyte transendothelial migration (Fig. 3b,c); these are composed of genes identified above (Fig. 3a). Methods such as the t-test do not have the sensitivity to capture these biologically relevant mechanisms through analysis of changes in gene expression.

Variance and RNA-seq methods

Analog microarray methods for studying gene expression are being replaced by RNA-seq methods¹⁵. However, such digital sequencing methods have a variance structure similar to that of analog microarray data. This concern is illustrated by an example in which the relationship between variance and average expression is analyzed for human kidney samples (two replicates) subjected to RNA-seq by Illumina’s next-generation sequencing method¹⁶ (Fig. 4a). The strong relationship between variance and mean expression demonstrates the failure of the normal assumptions that residual variances are normally distributed and hence are uncorrelated. Change presented as a function of expression in human liver versus kidney demonstrates the importance of analyzing the variance structure for the digital gene-expression data (Fig. 4b), similar to the results presented above (Fig. 3). The DESeq approach¹⁷, like VAMPIRE, exploits the expression-variance relationship of assay-wide data to better approximate gene-expression variances for experiments with few replicates. Similar to the microarray data sets presented above (Figs. 2a and 3a), transcripts with differences in expression are systematically distributed so that more confidence is given to measurements with a high sequence-tag count, whereas a greater change (higher ‘-fold’ value) is required for confident determination of the significance of features with a low sequencetag count.

Summary

Gene-expression measurements have become integral tools in modern immunology¹⁸. Given the paucity of samples, especially in human studies, it is essential that analysis methods help experimental designs with minimal sample size and number. The use of mean and fold-change ratios alone can be misleading and can lead to erroneous assignment of mechanistic changes or can miss important mechanisms that can help diagnostic and therapeutic processes. There is rich information in both analog and digital expression data, and use of the inherent error structure of data, as well as analysis of the variance and its use in modeling expression changes, can provide valuable information about markers and mechanisms of disease. Methods such as VAMPIRE, in addition to providing accurate analysis of samples of a limited size, can also be sufficiently sensitive to identify genes with different expression in regions of low tag count with low average expression. Most importantly, such methods avoid the need for arbitrary choices of fold-change cut-offs or assessment of significance with the t-test, which suffers from low statistical power.