alexa Screening for Differentially Expressed Genes: Are Multilevel Models Helpful? | Open Access Journals
ISSN: 2155-6180
Journal of Biometrics & Biostatistics
Like us on:
Make the best use of Scientific Research and information from our 700+ peer reviewed, Open Access Journals that operates with the help of 50,000+ Editorial Board Members and esteemed reviewers and 1000+ Scientific associations in Medical, Clinical, Pharmaceutical, Engineering, Technology and Management Fields.
Meet Inspiring Speakers and Experts at our 3000+ Global Conferenceseries Events with over 600+ Conferences, 1200+ Symposiums and 1200+ Workshops on
Medical, Pharma, Engineering, Science, Technology and Business

Screening for Differentially Expressed Genes: Are Multilevel Models Helpful?

Dongmei Liu1, Giovanni Parmigiani2* and Brian Caffo1

1Department of Biostatistics, Johns Hopkins School of Public Health, USA

2Department of Biostatistics and Computational Biology, Dana Farber Cancer Institute, Boston, USA

*Corresponding Author:
Giovanni Parmigiani
Department of Biostatistics and Computational Biology
Dana Farber Cancer Institute, 450 Brookline Avenue
Boston MA 02215, USA
Tel: (617) 632-5323
Fax: (617) 632-2444
E-mail: [email protected]

Received date: March 10 2014; Accepted date: April 09, 2014; Published date: April 15, 2014

Citation: Liu D, Parmigiani G, Caffo B (2014) Screening for Differentially Expressed Genes: Are Multilevel Models Helpful? J Biomet Biostat 5:192. doi: 10.4172/2155-6180.1000192

Copyright: © 2014 Liu D, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are are credited.

Visit for more related articles at Journal of Biometrics & Biostatistics

Abstract

Screening for changes in gene expression across biological conditions using high throughput technologies is now common in biology. In this paper we present a broad Bayesian multilevel framework for developing computationally fast shrinkage-based screening tools for this purpose. Our scheme makes it easy to adapt the choice of statistics to the goals of the analysis and to the genomic distributions of signal and noise. We empirically investigate the extent to which these shrinkage-based statistics improve performance, and the situations in which such improvements are larger. Our evaluation uses both extensive simulations and controlled biological experiments. The experimental data include a socalled spike-in experiment, in which the target biological signal is known, and a two-sample experiment, which illustrates the typical conditions in which the methods are applied. Our results emphasize two important practical concerns that are not receiving sufficient attention in applied work in this area. First, while shrinkage strategies based on multilevel models are able to improve selection performance, they require careful verification of the assumptions on the relationship between signal and noise. Incorrect specification of this relationship can negatively affect a selection procedure. Because this inter-gene relationship is generally identifiable in genomic experiments, we suggest a simple diagnostic plot to assist model checking. Secondly, no statistic performs optimally across two common categories of experimental goals: selecting genes with large changes, and selecting genes with reliably measured changes. Therefore, careful consideration of analysis goals is critical in the choice of the approach taken.

Background

Many genomics investigations using expression arrays take the form of searching for genes whose expression level is different across experimental conditions or phenotypes. The list of gene transcripts produced by a microarray analysis is usually the starting point for extensive additional biological work, including independent validation, and both in-silico and laboratory work on sequences and proteins related to the transcripts selected. In this context, microarray experiments are screening, not testing, experiments. Because of the wide range of important questions that can be explored using these arrays, and the costs involved, comparisons across conditions are often made using a limited number of replications. Efficient use of data is critical in improving a laboratory’s ability to correctly identify important biological hypotheses and proceed to test them by appropriate further experimentation.

Specific screening goals vary with the study. Two simple but representative situations are the selection of genes that are changed by a large amount, and the selection of genes that are changed by a reliably measured amount. In either case, the comparison of gene expression across two conditions based on replicated experiments requires a tradeoff of signal, the variation of expression across the two conditions, versus noise, the variation of expression within each condition. Therefore the problem is statistical in nature [1-3]. In this paper we discuss a Bayesian multilevel framework for developing screening tools that adapt to the goals of the analysis and to the genomic distributions of signal and noise. We evaluate a representative set of these tools using both extensive simulations and controlled biological experiments in which the set of altered genes is known.

A variety of approaches for selecting differentially expressed genes have been proposed [4] for a review and Murie et al. [5]. The simplest and still the most widely used is to set a threshold on a measure of signal alone, for example an estimated fold–change. This can be motivated by the desire to identify large changes, although often it is used by simple analogy with other gene expression essays that have much less noise. Upper and lower thresholds of two and one half are often seen in applications.

One limitation of this approach is that it does not consider how reliably gene-specific changes are measured. That is, it implicitly assumes that all genes are subject to the same level of noise. This may not be the case because even after appropriate preprocessing of the data, the within–gene variation in expression can be highly gene–dependent.

A straightforward way to account for both signal and noise is to select genes based on statistics motivated by two–sample testing, such as the T–ratio or the Wilcoxon statistic. For each gene, the T–ratio is an estimate of the signal–to–noise ratio. Because it requires estimating two or three parameters instead of one, when the number of replicates is small, the T–ratio does not necessarily perform better than fold– change, even when the goal is point-null-like. Gains in efficiency over both fold–change and T–ratios can be obtained by considering the ensemble of gene expression measures at once, rather than each gene in isolation. This occurs for at least two reasons. First, genes measured on the same array type in the same laboratory are all affected by a number of common sources of noise. Secondly, many changes in expression are part of common biological mechanisms.

A widely used approach that uses genome–wide information is “significance analysis of microarrays” or SAM [6]. SAM involves transforming the signal–to–noise ratios so that they are approximately independent of noise across genes. The type of transformation used by SAM is designed to protect against false discoveries generated by very small denominators in the T–ratios. The denominators that are really small are highly likely to be so by chance, because genes share many sources of variability, and a certain amount of variation is to be expected from all of them.

More broadly, joint estimation of many related quantities is often approached by multilevel modeling, and the associated Empirical Bayesian [7,8] and Hierarchical Bayesian [9] estimation techniques. See Carlin and Louis [10] for a detailed discussion. In genomics, these may represent variation in two stages. The first stage defines summaries at the gene level, for example test statistics, or estimates of fold change and noise. These describe variability of samples within each gene.

The second stage posits a “genomic” distribution for these gene–level summaries. Such multilevel modeling provides tools for borrowing strength from other genes when making inference on each gene. Some examples of implementations in microarrays are provided by Baldi and Long [11], Newton et al. [12], Efron et al. [13], L¨onnstedt and Speed et al. [14], Ibrahim et al. [15], Parmigiani et al. [16], Wright and Simon [17] among others.

In practice, a question often raised by genomics practitioners is the extent to which simple, real-time, shrinkage statistics motivated by multilevel models would outperform single-gene-at- a-time analysis or SAM. In this paper we set out to systematically address this question. To this end, we found it necessary to develop a general framework for developing and evaluating these fast shrinkage statistics. Our answer will turn out to be that shrinkage can furnish substantial improvement over single-gene-at-a-time analysis or SAM, provided that the statistics chosen will a) take into account the goals of the screening experiments and b) will be chosen based on examination of the properties of the genomic distribution of signal and noise. Even though we considered an extensive collection of statistics, the goal was not that of providing an exhaustive comparison of all approaches that have been proposed, but rather that of highlighting the critical role of the signal–to– noise trade–off and of providing tools to choose among alternative approaches based on the genome–wide behavior of signal and noise. We compared the performance of these statistics using both extensive simulations and real data sets in which fold changes were known.

Methods

Multilevel models for two-group comparisons

We consider a design in which two biological types are compared on a microarray that probes G genes. Each type is measured on n arrays using either technical or biological replicates. Here technical replicates refer to experiments that have multiple aliquots of the same RNA, while biological replicates refer to experiments that have multiple subjects from a population. Each situation requires a different interpretation of the array-to-array variability, but the formal structure is the same. We do not consider both levels of replication at the same time here. We denote by X1gj the expression for gene g in sample j in the first group, and by X2gj the expression for gene g in sample j in the second group. Expression levels are assumed to be centered around an overall experiment-wise mean.

Recall our interest lies in studying approaches for selecting genes that are differentially expressed between groups. We begin by describing an additive group effect and independent Gaussian errors. That is we assume that the observed expressions are conditionally independent draws from

equation

equation

Here, δg is the difference in expression level for gene g across groups, μg is an overall expression level for gene g, also referred to as abundance, or intensity, and equationis the variance of expression level for gene g in both groups. We refer to δg as true signal, and to σg as true noise. In a multilevel setting, our parameterization is different from that assuming E{X1gj}=μgE{X2gj}=μgg, which would lead to two different marginal variances in the two groups.

The fit of the normal distribution can often be improved by a suitable transformation of the data. Departures from normality, equal variances, and additivity may occur but are not considered in this manuscript.

Multilevel models postulate a distribution for the abundance, signal, and noise parameters across genes. A common assumption to many of the multilevel models used in microarray analysis is that of conjugate distributions for the second stage of the statistical model, in short a “conjugate model”. In the case of Gaussian data, the conjugate model implies that the gene–specific signal–to–noise ratios and abundanceto- noise ratios are independent of the corresponding gene–specific noise [18,19]. This assumption leads to convenient mathematical representations for many of the steps required by the data analysis, and is sometimes adopted solely for this reason. In practice, however, some microarray experiments follow this independence pattern closely, while others depart from it substantially. The loss of efficiency of screening based on the conjugate model in the latter case can be large.

Here we broaden the conjugate scheme and we investigate four model varieties, that result from the combination of two factors: (i) whether the gene-specific signal is independent of the gene-specific noise, (ii) whether the gene-specific abundance is independent of the gene-specific noise. Formally, for (i) the independence models assumes that δg and σg are independent, while the conjugate model assumes that δgg and σg are independent.

The remainder of our distributional assumptions are standard for normal multilevel models [9,20]. The models are summarized in Table 1.

Data Level

II. Independence
CI. Independence of Signal and Noise
IC. Independence of Abundance and Noise
CC. Complete Conjugacy

Table 1: The four classes of multilevel models investigated. The array-to-array variation is modeled in the same way in all four cases. In all cases, j = 1, 2, . . . , n and g = 1, 2, . . . ,G. All quantities denoted by Greek letters are unknown. A further set of prior distributions for the hyperparameters is described in the text.

We use the notation dg for the mean difference of expression across two groups, ag for the overall mean expression, and sg for the pooled estimate of the standard deviation. Notationally:

equation

equation

equation

In all four models of Table 1, these statistics are independent conditional on gene-specific parameters and have distributions

equation

equation

equation

conditional on gene-specific parameters. All four combinations of Table 1 occur commonly in practice. For example, our two experimental data sets show two markedly different relationships between signal and noise.

The vector of unknown genome-wide parameters will be denoted by ξ=(ν,β,λ,τ). There are several estimation approaches available for models of this kind. State-of-the art, computationally intensive approaches are usually based on MCMC [21]. Instead we focus on a faster and simpler empirical Bayesian approach based on estimating ξ by method of moments from the empirical distributions of equation, dg’s. The resulting estimators are computationally cheap and may include shrinkage of the signal, of the noise, or both. Several method of moments alternatives are available, and results can be strongly affected by this choice. For example, in our experience, the method of moments applied to the distribution of equation, which is inverse gamma, performs poorly, while the same applied to equation performs well.

We approach the task of generating a list of candidate genes by ranking genes according to a one-dimensional statistic, and then selecting all genes whose statistic is above a certain cutoff. This is the norm in practice. While more general decision theoretic approaches evaluating the trade-off between false and missed discoveries are available [22,23], these are complex, and would have been prohibitive in our extensive simulation study. The cutoff is often determined by the ability of a laboratory to perform validatory analyses, or, more inferentially, by false discovery rates [24-27]. In our presentation, to simplify the comparison of approaches, we focus on the ranking of genes implied by the statistics, and the ability of each statistic of identifying the top g genes.

Throughout, we draw a distinction between the selection of genes that are changed by a large amount, and genes that are changed by a reliably measured amount. Accordingly we consider two broad families of statistics, ones that estimate the signal, δg, and ones that estimate the signal-to-noise ratio, δgg. Because we use statistics as ranking devices and compare them based on ROC curves, we only need to define statistics up to constants that are not gene-specific. A proportionality sign will indicate omission of such constants. Table 2 summarizes the statistics we examined, organizing them by goal and motivating model structure. Table 3 summarizes the expressions of statistics motivated by multilevel models.

MOTIVATING MODEL ANALYSIS GOAL
Large change Reliably measured change
Independence/Normality of Genes Difference in Expression (F) T–statistic (T)
Exchangeability of Genes Significance Analysis of Microarray (SAM)
Complete Conjugacy Signal (CC.F) Standardized Signal (CC.T)
Tail probability (CC.TP) Bayes factor (CC.BF)
Independence of Abundance and Noise Signal (CI.F) Standardized Signal (CI.T)
Tail probability (CI.TP) Bayes factor (CI.BF)
Independence of Signal and Noise Signal (IC.F) Standardized Signal (IC.T)
Tail probability (IC.TP) Bayes factor (IC.BF)
Independence Signal (II.F) Standardized Signal (II.T)
Tail probability (II.TP) Bayes factor (II.BF)

Table 2: Summary of statistics examined, by goal and motivating model structure.

Model Statistics
F T BF TP
CC.m dg
CC.c dg
IC dg
CI
II

Table 3: Summary of functional forms of all the statistics motivated by multilevel models. For the complete conjugacy case (CC) we consider both statistics that use analytic integration with respect to σg (labelled CC.m) and statistics that use plug-in estimates of σg (labeled CC.c).

Statistics

In this subsection we enumerate and briefly comment on each of the statistics we considered. The remainder of this section is provided as a reference for future sections. Details of the derivations are given in the Appendix.

Difference in Expression (F): This is the observed average difference dg. Usually expression data are analyzed in the logarithmic scale, in which case F corresponds to an estimate of the log fold change across conditions.

T–statistic (T): This is the common statistics T ∝ dg/sg used for testing the null hypothesis of δg=0 one gene at the time.

Significance Analysis of Microarrays (SAM): This was proposed by Tusher et al. [6] and is based on the change of gene expression relative to an adjusted standard deviation. For the two group case considered here, the SAM statistic for gene g is

equation

where s0 is the so-called “exchangeability factor”. This factor is estimated using information from the entire set of genes to transform the values of SAM so that noise and SAM are approximately independent.

Statistics for the Complete Conjugacy (CC) Model: In the Complete Conjugacy model the conditional posterior distribution of δg given hyperparameters ξ can be written as

            equation                       (1)

A set of computationally cheap statistics is derived by considering empirical Bayes estimates of this distribution, obtained by replacing hyperparameters ξ with an estimate equation. A regularized estimate of signal δg is

equation

Regularization is independent of the gene, so for any given experiment CC.F will be proportional to F. For this reason we only consider F, although we keep this correspondence in mind when interpreting the results.

A standardized estimate of signal is derived as the ratio of the conditional posterior mean and standard deviation of δg from expression (1),

equation

The denominator incorporates a linear shrinkage estimate of the gene-specific variance with gene-varying coefficients, penalizing more heavily genes whose signal or abundance are outlying. For this reason, it is critical that the conjugacy assumption be checked, or very valuable information may be lost. On the other hand, when the assumption is met, an increase in efficiency is gained from estimating the denominator.

An empirical Bayes estimate of the Bayes factor [28] for the null hypothesis of no gene-specific differential expression, δg = 0, is:

equation

Finally, we consider the empirical Bayes approximation, equation of the probability that the true change δg exceeds D [29]. Here, D represents a target change across conditions. This tail probability reflects the observed change, its variability and the likely magnitude of biologically significant changes.

Statistics for the Independence of Abundance and Noise (IC) Model: In this model the posterior distribution of δg given σg and hyperparameters ξ can be written as

            equation                       (2)

Unlike in the complete conjugate case, a closed form marginalization with respect to σ is not possible. Therefore we derive results assuming equationis known. In the actual calculations, to obtain a real-time statistics, equationis estimated by the posterior mode of the distribution of equation , ξ. Then our estimate of the normalized signal is

equation

As with CC.F, IC.F is proportional to F, so we only consider F in our results section.

A standardized estimate of signal based on regularized estimates of signal is the ratio of the marginal posterior mean and standard deviation of δg from expression (2), that is

equation

The Empirical Bayes estimate of the Bayes factor, conditional on gene specific variance is:

equation

Finally, we consider the empirical Bayes tail probability

equation

Statistics for the Independence of Signal and Noise (CI) Model: In this model the posterior distribution of δg given σg and hyperparameters ξ can be written as

            equation                       (3)

Again, to obtain a real-time statistic we develop results conditional on equation and estimate it with its posterior mode in actual calculation. The estimate of the signal is

equation

A standardized estimate of signal based on regularized estimates of signal is the ratio of the marginal posterior mean and standard deviation of δg from expression (3), that is

equation

The Empirical Bayes estimate of the Bayes factor is:

equation

while the tail probability approximation is

equation

Statistics for the Complete Independence (II) Model: In this model the posterior distribution of δg given σg and hyperparameters ξ can be written again as

            equation                       (4)

A regularized estimate of δg, motivated by the independence model is obtained by replacing ξ with equationand equationwith its conditional posterior mode evaluated at equation , and approximating the posterior mean by

equation

Unlike IC.F and CC.F, both II.F and CI.F imply a linear shrinkage which depends on the genomic variability of the signal. Dividing II.F by the square root of the variance of δg, and approximating as before, we obtain

equation

The Empirical Bayes estimate of the Bayes factor is:

equation

while the empirical Bayes approximation to the tail probability is

equation

Notice that the definitions of statistics for the CI and II cases would be the same if equation was known. Hence the only difference in practice is the posterior mode for equation. . It should not then be a surprise that the performance of these two are very close. For the same reason, this is also true for statistics in the CC and IC settings conditional on equation.

The empirical Bayes estimators, both standardized and not, have functional similarities to the SAM score, although shrinkage of the noise in the denominators are determined differently. In empirical Bayes analyses, the shrinkage is driven by the parameters of the genomic distributions of signal and noise, in a form that depends on whether or not conjugacy is assumed. In SAM one applies linear shrinkage to the standard deviation rather than the variance, and the shrinkage intercept s0 is chosen to approximate independence of SAM ratios from noise.

Simulation Results

Overview

In this section we study the performance of the real time shrinkage statistics on a large number of data sets simulated from each of the four models in Table 1. We evaluate each statistic on the basis of the implied ranking of genes, and the ability of each statistic of identifying the top g genes. In the analyses presented here we considered two alternative goals: in one the genes of interest are the top genes by absolute change δg. In the other the genes of interest are the top genes by signal-to noise ratio δgg. We considered the top 1%, 2% and 10% for each goal. For each cutoff we create a binary indicator of whether the true parameter is in the top list, and use this indicator as the true class assignment to be predicted. We evaluate performance by an ROC curve [13,30], which is the graph of the true positive fraction versus the false positive fraction for varying thresholds.

As summaries, we consider the overall area under the ROC curve [31] and the partial area under the ROC corresponding to false positive fractions smaller than 20%. We prefer these measures to others incorporating explicitly δg and δgg for two reasons: the goal of the microarray experiments we are focusing on is screening rather than estimation; interest usually lies in a relatively small fraction of important findings.

Based on these criteria, our simulations suggest three general conclusions about the alternative approaches for identifying differential genes: i) simple, real-time, shrinkage statistics motivated by multilevel models can outperform alternatives based on analyzing each gene separately, in some cases by a large margin; ii) the same statistics can perform better than the commonly used SAM [6] statistic, provided that careful checking of the multilevel modeling assumptions is carried out, and iii) no statistics is optimal for both the identification of large δg and large δgg.

Design of simulation study

Our goal was to generate a large and diverse number of scenarios, depending on sample size, conjugagy assumptions, and hyperparameter choices. We considered three sample sizes: 3, 10 and 100 per group. Use of samples as small as 3 is a common scenario in the gene screening experiments taking place during the routine activities of many laboratories, while 10 per group is a common scenario in comparisons across conditions for population genomic studies. Sample size as large as 100 per group are rare and considered here mostly as a check.

For each combination of conjugacy assumption (CC, IC, CI and II) and sample size, we simulated data from 2009 hyperparameter combinations, resulting in a total of 24108 datasets. The 2009 combinations of hyperparameters are based on the grid:

equation

equation

equation

equation

For each combination we derive ν and β from equationand equation. Here the total number of combinations is 2009 rather than the full 2401 because some expectation/variance combinations lead to unrealistic settings for _ yielding numerically unstable results.

Simulation results

Mining the massive information generated by the thousand datasets required drastic summarization. Here, we present one detaset in detail, and then provide the following summaries: scatterplots in which each deatset/statistic combination is represented by a single point, summaries of pairwise comparisons of statistics by model/goal/ sample size, and summaries of best performing statistics by model/ goal/sample size.

Figure 1 shows the ROC curves for a single simulation. Data are generated from the II model and the ROC is based on identifying genes with large signal. In this data set, we see a clear separation in the performance of the statistics, both within and across conjugacy structures. The tail probability statistics perform best irrespective of the motivating model, stressing the importance of correctly specifying the analysis goal.

biometrics-biostatistics-roc-curves

Figure 1: ROC curves for all the statistics under study in a single simulation with 3 replicates and 1000 genes. The hyperparameter values used to simulate the data were ν=2, β=1, γ=1 and τ=1. The ROC curves are grouped by statistical model to unclutter the displays. The curves for T and SAM are repeatedin each panel.

To summarize this type of comparison for all 2009 simulated data sets, we display pair scatterplots of areas under the ROC curves (AUC). Figure 2 summarizes results for II data with three replicates. We focus on T, SAM and on the three shrinkage-based statistics that perform best in II data. Whiskers at the top of the graphs for the T statistics indicate that, for a subset of simulation for which other statistics achieve a perfect separation, the T is can still miss a fraction of differentiated genes. The reciprocal situation does not occur, suggesting a substantial inefficiency in the use of the T statistic. For both goals, the best performing model is based on CI, suggesting that the model chosen to represent the abundance/noise relationship is less critical than that chosen for the signal/noise relationship. A complete set of such comparisons including IC, CI and CC is available in the supplementary materials on our website.

biometrics-biostatistics-comparison-of-performance

Figure 2: Comparison of performance of T, SAM and the shrinkage statistics that perform best in independent (II) data. Each point represents the areas under the curve for two statistics for a particular simulated data set. In the plots above the diagonal, we report the areas when identifying genes with high signal to noise ratio, while in plots below the diagonal we report areas when identifying genes with large signal. The shrinkage statistics that perform best change across the two goals, and therefore different shrinkage statistics are shown above and below the diagonal. Because over-plotting points may lead to visually misleading results in some cases, we also print the percentages PU and PL of points lying above and below the diagonal, and the averages MX and MY of the horizontal and vertical variables.

Rather than reproducing these plots for each of the simulation scenarios, we present, in Figure 3, “heat maps” synthesizing pairwise comparisons of estimators for the four simulation settings. For each pair of statistics, maps encode how often one statistics’ AUC is better than the other’s, allowing for rapid comparison of any two statistics on data generated under each of the four model families. To further summarize, the best performing statistics from the heat maps for sample sizes 3, 10 and 100 replicates are reported in Tables 4 and 5 (for partial AUC).

biometrics-biostatistics-heat-maps

Figure 3: Heat maps of the results for data simulated from the four models, with sample size 3. The color encodes the sign of the difference between percentages PU and PL of points lying above and below the diagonal in Figure 2 for all pairs of statistics. Green indicates PU<PL, red indicates PU>PL, while yellow indicates those cases where the difference was less than 0.05. The best statistics to identify genes with large signal are labeled with purple, while the best statistics to identify reliably measured differentially expressed genes are labeled with blue. These results are further summarized in Table 4.

biometrics-biostatistics-best-performing-statistics

Table 4: Summary of best performing statistics by AUC. Here “model” is the true model used for simulation. Statistics were evaluated by their ability to identify the top 1%, 2% and 10% of genes with large signals (labeled Signal) and large signal-to-noise ratios (S/N). The rows labeled MM correspond to parameter esti-mation using the method of moments while those labeled TT correspond to using the true hyperparameter. Instances where the best performing statistic did not match the appropriate true model and goal are high- lighted in red. Instances where the best performing statistic did not match the model but was consistent with goal are highlighted in blue. A wildcard * indicates the statistic from either the conjugate or indepen- dence model was the best performer. For example, ”*C.F” indicates that the CC.F and IC.F statistics were roughly equivalent best performers.

biometrics-biostatistics-best-statistics-simulation

Table 5: Best statistics for each simulation scenario with three replicates. Here “model” corresponds to the true model used for simulation. Statistics were differentiated in their ability to identify the top 1%, 2% and 10% of genes with large signals (labeled Signal) and large signal-to-noise ratios (S/N). The rows labeled “MM” correspond to parameter estimation using the method of moments. Results using the actual true parameter values (labeled "Truth”) are also given. Instances where the best statistic did not matching the appropriate true model and goal are highlighted in red. Instances where the best statistic did not match the model but was consistent with goal are highlighted in blue. The results highlight the importance of matching the statistic to the goal; in only one case did a T statistic out-perform other statistics to identify high signal changes while in no cases did a fold change out-perform other statistics to identify high SNRs.

Overall, these results emphasize the importance of matching the statistic to the analysis goal: fold changes and tail probabilities appear to be the best statistics for estimating large signals changes. In contrast the signal-to-noise ratios and Bayes factors appear to be optimal for estimating reliably measured differential expression. Also, these result confirm the importance of shrinkage: in only one case did a T statistic out-perform other statistics to identify high signal changes while in no cases did fold change out-perform other statistics to identify high SNRs. IN shrinkage correctly identifying the appropriate modelling assumptions becomes increasingly important as the number of replicates increases. The SAM statistic performs well across models, though often worse than the best shrinkage statistics within a model.

To provide a bound on the improvement in performance that can be achieved by shrinkage, in Tables 4 and 5 we provide results obtained by plugging in the true values of the parameters of the genomic distributions instead of their estimates in the calculation of the statistics.

When there is no close form, we have used the conditional estimation by plugging in the posterior mode of equation. To verify how reasonable these real-time statistics are, we considered the Complete Conjugacy model where closed forms are available for all the statistics. We performed another set of simulation with exactly the same hyperparameters comparing the results based on CC.TP and CC.BF in two cases: one based on the conditional posterior of δg with posterior mode of equation plugged in, the other based on the marginal posterior of δg integrating out equation. Results are shown in the supplementary materials. Real-time statistics generally perform as well as the exact statistics, although there is a minority of cases in which the exact statistics does far better than real-time approximation, especially for tail probabilities.

Finally, we further investigated the seemingly counterintuitive result where the best performing statistic for the data simulated from the Complete Conjugacy model are the IC.TP, IC.BF for genes ranked by both signal alone and signal-to-noise ratio. This behavior persists at larger sample sizes. The reason for the counterintuitive behavior is that some of the hyperparameter combinations lead to simulated dataset that have diagnostic plots consistent with an IC model, in which case IC statistic performs well while the abundance-based shrinkage applied by the CC statistics leads to loss of some of the signal. Additional details are provided in the supplementary materials.

Experimental Results

Datasets

We analyzed two data sets. The first was reported by Tusher et al. [6] in the context of comparing radiated and unirradiated cell lines. A subset of the genes’ changes, identified based on the SAM statistic, were subsequently validated by independent essays. While the experiment includes some blocking, we analyze it here as though it were a two– class comparison with 4 replicates.

The second data set is from an experiment reported by Dudley et al. [32]. They performed a so called “spike-in” experiment in which they selected 9 genes with very low natural expression and “spiked-in” Cy3- labeled gene-specific oligonucleotides in increments from 0.5 fold to 200 fold. Their experiment used cDNA microarrays including a total of 6307 genes, and had two replicates. We work from ratios of Cy3- to-Cy5 channels, after normalization [33]. While spike in experiments are useful in that true fold changes are known, both the magnitudes of the changes, and the sparsity of changes in the genome are unlikely to be realistic.

Graphical diagnostics for conjugacy structure

We begin by investigating the relationship between signal, abundance and noise by displaying boxplots of signal and abundance by noise level. These elaborate on ideas of Dudoit et al. [34] and Tusher et al. [6]. Figure 4 considers the Tusher data. An SN plot in which the location and dispersion of signal are stable across noise levels suggests the use of an independence model, while one in which the location and dispersion of signal increase with noise level suggests the use of a conjugate model. Thus, a constant box size indicates independence while an increasing box size indicates conjugacy. In simulated data (see Supplementary materials), the diagnostic plot clearly distinguishes conjugacy with respect to signal and noise as well as conjugacy with respect to abundance and noise.

biometrics-biostatistics-signal-to-noise-diagnostic

Figure 4: Signal-to-noise diagnostic plot for the Tusher data under two transformations. Genes identified by SAM are highlighted in red. The horizontal axis is the binned gene specific variances. The vertical axis for the plots in the left column is the average expression across the two groups (abundance) while it is the the average difference in expression (signal) for the right column. Conjugacy appears appropriate for the raw data, and an independence relationship appears appropriate for the transformed data.

Because results are sensitive to the type of transformation applied to the expression measurements, we display both the original scale and the cube root. Untransformed data show a pattern consistent with the conjugate model, while data transformed using the cube root appears consistent with the independent model. An alternative visualization to the SN plot is a simple scatterplot of signal versus noise. A limitation of this approach is that it can be difficult to establish whether increased variation in signal at different level of noise is due to a true relationship or simply to a higher number of genes at that noise level.

Figure 5 shows the spike-in data using three transformations. The original scale shows a marked positive relationship between estimated signal and noise, the cube root scale a mild positive relationship, and the logarithm an almost stable relationship, with some indication of larger variation in the signal at lower noise level. These figures do not inform us about absolute intensity, so the larger variation of signal at the low end after the log transformation is not the same as the wellknown “fishtail” effect observed in MVA plots.

biometrics-biostatistics-signal-to-noise-plot

Figure 5: Signal to noise plot for the spike-in data under three transformations. Spiked-in genes are high- lighted in red. The horizontal axis is the binned gene specific variances while the vertical axis for the plots in the left column is the average expression across the two groups (abundance) while it is the the average dif- ference in expression (signal) for the right column. The apparent relationships between abundance and noise and signal and noise clearly change dependent on the transformation used. While conjugacy appears ap-propriate for the raw data, and independence relationship appears more appropriate for the log-transformed data.

In this data set the cube root transformation is the most effective in helping identify the truly differentially expressed genes, performing better than the commonly used log. In practical application one does not have the advantage of knowing the true changes when choosing a transformation. The important lesson here is, however, that choosing transformations based on convenient statistical properties such as variance stabilization does not necessarily improve, and could prejudice, our ability to detect signal.

These two datasets stress that both the independence of signal and noise and independence of signal-to-noise ratio and noise may need to be tackled in real applications. While transformation of the measured intensities may allow one to achieve independence, it is not clear that such transformations would be optimal in terms of gene screening.

Comparison of real-time statistics

Figure 6 shows pairwise scatterplots of the statistics CC.TP, II.F, and SAM for the two transformations in Figure 4. CC.TP and II.F are the two best performing real-time shrinkage statistics for identifying genes with large signal. The two best statistics for selecting reliably measured genes, CC.BF and II.T, are given in the supplementary materials. All the genes originally identified by Tusher et al. [6] receive high tail probability using both transformed and untransformed data, though additional genes also receive tail probability close to one. On the other hand, the correspondence between SAM and II.F is good after cube root transformation but not in the raw scale. In evaluating these results, one must keep in mind that only genes that exceeded a certain SAM threshold were validated independently in the study. Therefore, direct performance comparisons with SAM are not reliable here.

biometrics-biostatistics-tusher-data

Figure 6: Comparisons of the II.F, CC.TP and SAM statistics on the Tusher data. Genes selected by SAM are highlighted in red for reference. Here the CC.F and II.TP statistics were identified from Tables 4, 5 and 6 as the optimal statistics for detecting large signal changes for the II and CC models. For this data set, the II model is supported under the cubic transformations while the CC model is supported for the raw data.

Figure 7 compare statistics in the spike-in data. For the log transformed data, we would expect a better performance from II.F than CC.TP based on the SN plot. In fact, the II.F statistics shrinks the effects excessively and gives a less efficient ranking. For the cube root transformed spike-in data we would expect and, in fact, see a better performance from CC.TP than II.F in Figure 7. For the untransformed spike–in data, results, shown in the supplementary materials, confirm the intuition from the exploratory plots that the conjugacy model should outperform the independence model. Figure 7 does show this result. In a close view of comparison between CC.TP and SAM on raw data (see Supplementary materials), CC.TP also clearly picks up all the spiked genes, while SAM does not. Spiked genes are genes have large signals, so that poor performance of II.T and CC.BF is no surprise (see Supplementary materials).

biometrics-biostatistics-fold-change-statistics

Figure 7: Comparisons of fold change statistics II.F and CC.TP with SAM in the spike-in data. Spiked-in genes are highlighted in red. Here the CC.F and II.TP statistics were chosen from Tables 4, 5 and 6 as the optimal statistics for detecting large signal changes for the II and CC models. For this data, the II model is supported under the log transformations while the CC model is supported for the raw and cubic transformed data. The II.F statistic shrinks conservatively even for the log transformed data while SAM and CC.TP perform well for both transformed data sets. No statistic performs well on the original scale.

This analysis warns of a potential danger in the use of multilevel modeling. Assuming that the signals follow a Gaussian law assumes that all genes are differentially expressed to some extent, and the goal is to either detect the largest signals or the largest reliably measured signals. This is realistic in case-control comparisons and in experiments in which the experimental intervention changes a large portion of the expression, as during cell division. In contrast, in this spike-in experiment, the distribution of signals is in fact degenerate at 0 for all but the 8 spiked in genes. Therefore, the statistics motivated by the Gaussian law on the signal are not validated from the data. While extreme, the spike-in situation may be relevant in practice when experimental intervention modifies a small set of genes involved in a very specialized pathway.

Diagnosing empirically whether the signal distribution is a mixture is difficult. Appropriate weight should be given to the biological circumstances of the experiment. For example, here an independence relationship is suggested for the signal and noise for the cubic and log transformed data. However, as the majority of the genes are biologically known to have no signal, these plots do not inform us on the question of interest. Furthermore, Figure 7 shows that the best complete independence statistics for identifying large signal and reliably measured signal changes, II.F and II.T, perform poorly for detecting the spiked-in genes. In summary, aggressively modelling the distribution of signals when the overwhelming majority of genes have no signal can produce poor results.

Conclusions

In this article we present a framework for interpreting, selecting, and estimating shrinkage based screening statistics used in the identification of differentially expressed genes. We also evaluated a representative set of these tools using both extensive simulations and controlled biological experiments in which the set of altered genes is known or partially known.

Our results emphasize two important practical concerns that are not receiving sufficient attention in applied work in this area. First, while shrinkage strategies based on multilevel models are able to improve selection performance, they require careful verification of the assumptions on the relationship between signal and noise. Incorrect specification of this relationship can negatively affect a selection procedure. Because this inter-gene relationship is generally identified in genomic experiments, we suggest a simple diagnostic plot to assist model checking. Secondly, no statistic performs optimally across two common categories of experimental goals: selecting genes with large changes, and selecting genes with reliably measured changes. Therefore, careful consideration of analysis goals is critical in the choice of the approach taken.

The commonly used SAM statistics emerges as a reasonable compromise between the two goals above and is, to some extent, automatically adaptive to different relationships between signal and noise. Improving on SAM is possible but requires careful validation of the assumptions about the upper level distribution. The assumption of conjugacy in the abundance dimension requires careful attention as it is not robust. In particular, estimators based on the CC assumption can be outperformed even on data generated under CC.

Our simulation analysis relies on the assumed normality of data. In practice, two aspects of it are critical. At the lower stage, in small samples, the functional form of the error distribution across samples is hard to assess. On the other hand, at the upper stage,normality can be checked, and transformations may help, although the caveats discussed in Section 4.3 should be considered. Alternative multilevel models have been studied, for example by Newton and Kendziorski [35] who consider gamma models, and M¨uller et al. [22], who extend those to mixtures of gamma models. While these alternatives are worth serious consideration, here we focus on Gaussian models and statistics motivated by the Gaussian setting, primarily because in this way we can practically investigate a variety of relevant statistics on a massive number of simulation scenarios.

Other interesting multilevel approaches have been proposed to analyze designs that are more complex than the two-group comparisons considered here. We refer the reader to Kerr et al. [36], Wolfinger et al. [37], Kooperberg et al. [38], Tai and Speed [9] and Wang et al. [40] for further details. Meta-analysis of multiple microarray studies is another area where multilevel models have proven helpful, as illustrated by Conlon [41] and Scharpf et al. [42].

In recent years gene expression is increasingly measured using technologies based on sequencing short read (RNA-seq) which coexist with the hybridization-based approaches that motivated this work. These technologies generate count data. Though the normal model is sometimes used and can perform reasonably after variance stabilization, differential expression is better analyzed statistically using poisson and negative binomial models as described in [43] or Rapaport et al. [44]. An effective effective approaches using multilevel models is included in the Bioconductor package described by Love et al. [45].

In our analysis we assumed that all genes on the array are potentially changed. This is realistic in case control designs across populations or comparison of cells at different stages of the cell cycle, regulation can be expected in a large number of genes, although differences will vary randomly and many genes will be changed by amounts that are smaller than noise. In tightly controlled experiments, such as a comparison of wildtype versus mutant species, or treated versus untreated cell lines, differential expression may only involve a small number of pathways and genes. While this situation can be reasonably handled in the framework considered here, it would be more accurately modeled by assuming that only a fraction of the genes are differentially expressed across groups A multilevel model could assume assume that a fraction of the δg are identically zero, while the rest are normally distributed [14].

Our focus here has been on simple and easy-to-compute statistics for gene selection. Multilevel models were used only to provide motivation and a conceptual framework for the derivation of the shrinkage statistics. More generally, multilevel models give raise to potentially more efficient strategies than those considered here, at the price of increased computational expense, for example by MCMC of MCEM. Systematic exploration of those in the thousands of data sets considered here would have been impractical. However, our results suggest that that appropriate shrinkage is a critical part of gene selection and this will hopefully encourage practitioners to consider these more computing intensive approaches as well.

Acknowledgement

Data from Dudley, Aach, Steffen and Church (2002) was downloaded from http://arep.med.harvard.edu/masliner/supplement.htm. Data from Tusher, Tibshirani and Chu (2001) was downloaded from http://www-stat-class.stanford.edu/SAM/SAMServlet. Work of Giovanni Parmigiani supported by grant NCI 5P30CA006516-49.

References

Select your language of interest to view the total content in your interested language
Post your comment

Share This Article

Relevant Topics

Recommended Conferences

Article Usage

  • Total views: 11620
  • [From(publication date):
    April-2014 - Oct 17, 2017]
  • Breakdown by view type
  • HTML page views : 7855
  • PDF downloads :3765
 

Post your comment

captcha   Reload  Can't read the image? click here to refresh

Peer Reviewed Journals
 
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals
International Conferences 2017-18
 
Meet Inspiring Speakers and Experts at our 3000+ Global Annual Meetings

Contact Us

Agri, Food, Aqua and Veterinary Science Journals

Dr. Krish

[email protected]

1-702-714-7001 Extn: 9040

Clinical and Biochemistry Journals

Datta A

[email protected]

1-702-714-7001Extn: 9037

Business & Management Journals

Ronald

[email protected]

1-702-714-7001Extn: 9042

Chemical Engineering and Chemistry Journals

Gabriel Shaw

[email protected]

1-702-714-7001 Extn: 9040

Earth & Environmental Sciences

Katie Wilson

[email protected]

1-702-714-7001Extn: 9042

Engineering Journals

James Franklin

[email protected]

1-702-714-7001Extn: 9042

General Science and Health care Journals

Andrea Jason

[email protected]

1-702-714-7001Extn: 9043

Genetics and Molecular Biology Journals

Anna Melissa

[email protected]

1-702-714-7001 Extn: 9006

Immunology & Microbiology Journals

David Gorantl

[email protected]

1-702-714-7001Extn: 9014

Informatics Journals

Stephanie Skinner

[email protected]

1-702-714-7001Extn: 9039

Material Sciences Journals

Rachle Green

[email protected]

1-702-714-7001Extn: 9039

Mathematics and Physics Journals

Jim Willison

[email protected]

1-702-714-7001 Extn: 9042

Medical Journals

Nimmi Anna

[email protected]

1-702-714-7001 Extn: 9038

Neuroscience & Psychology Journals

Nathan T

[email protected]

1-702-714-7001Extn: 9041

Pharmaceutical Sciences Journals

John Behannon

[email protected]

1-702-714-7001Extn: 9007

Social & Political Science Journals

Steve Harry

[email protected]

1-702-714-7001 Extn: 9042

 
© 2008-2017 OMICS International - Open Access Publisher. Best viewed in Mozilla Firefox | Google Chrome | Above IE 7.0 version
adwords