A full Bayesian hierarchical mixture model for the variance of gene differential expression
© Manda et al. 2007
Received: 11 October 2006
Accepted: 17 April 2007
Published: 17 April 2007
Skip to main content
© Manda et al. 2007
Received: 11 October 2006
Accepted: 17 April 2007
Published: 17 April 2007
In many laboratory-based high throughput microarray experiments, there are very few replicates of gene expression levels. Thus, estimates of gene variances are inaccurate. Visual inspection of graphical summaries of these data usually reveals that heteroscedasticity is present, and the standard approach to address this is to take a log2 transformation. In such circumstances, it is then common to assume that gene variability is constant when an analysis of these data is undertaken. However, this is perhaps too stringent an assumption. More careful inspection reveals that the simple log2 transformation does not remove the problem of heteroscedasticity. An alternative strategy is to assume independent gene-specific variances; although again this is problematic as variance estimates based on few replications are highly unstable. More meaningful and reliable comparisons of gene expression might be achieved, for different conditions or different tissue samples, where the test statistics are based on accurate estimates of gene variability; a crucial step in the identification of differentially expressed genes.
We propose a Bayesian mixture model, which classifies genes according to similarity in their variance. The result is that genes in the same latent class share the similar variance, estimated from a larger number of replicates than purely those per gene, i.e. the total of all replicates of all genes in the same latent class. An example dataset, consisting of 9216 genes with four replicates per condition, resulted in four latent classes based on their similarity of the variance.
The mixture variance model provides a realistic and flexible estimate for the variance of gene expression data under limited replicates. We believe that in using the latent class variances, estimated from a larger number of genes in each derived latent group, the p-values obtained are more robust than either using a constant gene or gene-specific variance estimate.
The recent advancement of deoxyribonucleic acid (DNA) microarray technology allows the measurement of expression levels of tens of thousands of genes simultaneously [1, 2]. A DNA microarray experiment measures the abundance of messenger ribonucleic acid (mRNA) present in a set of cells, and a high concentration of mRNA for a given gene indicates a high expression level for that gene [3, 4]. The solution of mRNA is either radiolabelled or fluorescently labelled and then allowed to hybridize to spots on the array. Further information about early DNA microarray experiments can be found in . There is a wide variety of arrays, but the two main kinds are short and long oligonucleotide arrays. In oligonucleotide arrays, there are approximately 20 probes pairs of a perfect match (PM) and a mismatch (MM) for each gene. The PM probe contains a match to a small subsequence of a gene's polynucleotide, about 25 bases long, while the MM probe acts as a control, being a copy of the PM probe with the central position flipped to its complement. The amount of mRNA present in the target gene in a sample is estimated from a combination (usually, the average) of the PM-MM intensity differences over the 20 probe pairs .
In a spotted complimentary DNA (cDNA) or long oligonucleotide arrays (the kind featured in this article) thousands of spots of cDNA from the genes are printed onto a glass slide or some other form of substrate. Two different mRNA samples are separately reverse-transcribed into cDNA and labelled with different fluorescent dyes, green (cyanine 3) or red (cyanine 5). The mixture of labelled cDNA is co-hybridized on the same microarray, and the labelled cDNA molecule will bind to the complementary fragments of cDNA sequence on the slide. A laser scanner is then used to measure both fluorescent signals emitted at each spot on the chip. The general idea assumed behind the technology is that if a particular gene is highly expressed in the sample, it produces many molecules of mRNA . These in turn will hybridize to the probes on the microarray and generate a very bright fluorescent area. Genes that are less expressed produce less mRNA, which results in dimmer fluorescent spots. If there is no fluorescence, no messenger molecules have hybridized to the probes, indicating that the gene is inactive. By comparing the intensity levels of the emitted fluorescent lights between the samples, it is hoped that one might be able to identify any differences in the gene expression profiles of the various samples.
In a typical cDNA microarray experiment, we are looking to ascertain whether gene i displays differential expression between two samples T and C, labelled with differential colours red and green. For instance, these experiments include comparing tumor and normal tissue cells, treated and untreated cells, or cells from knockout and wild-type organisms. The samples can be compared on the same slide (i.e. same hybridization), resulting in a direct estimate of differences in expression levels since the measurements come from the same hybridization. An alternative is when expression levels T i and C i are estimated in two different hybridizations, with T together with reference sample R and C with another reference sample R'. This is an indirect estimate of the gene's differential expression since the T and C expression levels are from different hybridizations . Sometimes, common reference samples are hybridized with each mRNA sample of interest (T or C), resulting in what are known as common reference designs. The common reference sample could be tissues from wild-type organisms or control tissue or a pool of all the samples of interest.
The simplest way, used by many in the field, to ascertain a gene's differential expression, is on the basis of a fold-change criterion, defined by the log-ratio δ i = log2(T i /C i ), under direct comparison or δ i = log2(T i /R i ) - log2(C i / ) for indirect comparison. It is expected that the majority of genes will have a δ i value close to 0 . Those genes with a large positive δ i value (δ i > 1) are generally concluded to be overexpressed or upregulated in the T sample, and those whose δ i is negatively large (δ i < -1) are concluded to be underexpressed or downregulated in the T sample. However, the use of fold change is of limited use, as the intensities are associated with some biological, experimental and measurement error, and unless these error distributions can be derived, it is difficult to assess whether a ratio of 1.9, say, is worth noting or whether it has occurred by chance. Furthermore, the boundaries accepted for thresholding these fold-changes seem to be rather arbitrary and very little documentation can be found to support these criteria.
In recent times, the identification of genes that are differentially expressed between two conditions has been based on a standardised fold-change, which is the fold-change divided by an estimate of its standard deviation. A t-test, with a correction for multiple testing, is then used to test for significance of the standardised fold-change. This allows for the assessment of significance of the observed differences in the presence of all the sources of variation, which are not necessarily equal from gene to gene. Our contribution to the problem of identifying genuinely differentially expressed genes is on the estimation of reliable standard deviations of gene expression levels. Modelling of gene expression variability ranges between two extreme cases: a constant variance, which is too stringent an assumption, to independent gene-specific variances. The latter option has low power as it is based on very few replications as a result of the relatively large cost of commercial microarray chips. This makes estimation of the sample standard deviation very unreliable and unstable. An ad-hoc solution to the problem includes discarding genes with a small fold-change and very small standard deviations . A better method, called Significance Analysis of Microarrays (SAM), developed in , adds a constant a 0 to each gene-specific standard deviation, thus preventing the denominator of the t-statistic from getting too small. This was expounded in , using an empirical Bayesian procedure, by taking a 0 equal to the 90th percentile of the standard deviations of all the genes. Shrinking the gene-specific standard deviations in this way helps to minimize the false discovery rates (i.e. a large t statistic).
An intuitive approach to modelling gene expression data is to assume two groups of genes, one group with genes that are differentially expressed and the other with genes that are not differentially expressed. This approach has been used in analyses involving mixture models for gene expression levels. Formally, the basic assumption is that the distribution of the difference δ i can be flexibly modelled as a mixture with two components: a subgroup of genes with δ i around 0 and a subgroup of genes with non zero δ i [3, 4, 9, 10]. Using this approach, Lonnstedt and Speed , derived an empirical Bayesian statistic B, which is the log posterior odds of differential expression. Efron et al.  and Efron and Tibshirani  also consider a two component mixture model to model differential gene expression, the later using a rank-based nonparametric two-sample test statistic. A similar approach was followed-up in  based on a fully Bayesian hierarchical model, but with an unknown number of mixture components. The number of components was treated as a random variable and estimated with the other parameters based on the pioneering work of Richardson and Green . Other approaches using mixture models for gene expression data from microarray experiments can be found in [13, 14].
A number of methods, particularly based on full Bayesian hierarchical models, have been used to provide better estimates of variance for gene differential data. These methods provide estimates of gene-specific variance, which are the weighted average of the empirical variance and a prior variance estimate [11, 15, 16]. In particular, Lewin et al.  provide a fully Bayesian approach combining estimation of gene differential expression, biological variability and array effects with a hierarchical prior distribution on gene-specific variances. Other than modelling the gene-specific variances with an exchangeable hierarchical prior, Delmar et al.  use a mixture model on the distribution of gene-specific variances. Genes are grouped into latent classes based on homogeneity of their variances. A gene is assigned a variance based on its latent class membership and this variance is estimated with high power because of the large number of genes (hence a larger number of replicates) in that latent class. All these methods produce what are called regularized t-tests.
Our work is closely related to that proposed by Delmar et al. , who used the EM-algorithm to fit a variance mixture model for gene expression data. We believe our Bayesian approach has certain advantages and adds value in comparison to the EM-algorithm approach. In using a Bayesian hierarchical model, we are able to model various sources of variability in a common model, thus propagating uncertainty. Within a Bayesian hierarchical model framework, it is far easier to borrow and share data across all genes in order to obtain more reliable estimates of their variance and at the same time allowing for some variability. In this approach, variances are stabilized and shrunk towards the average variance within each component of the mixture, in particular some small and large variance estimates that are incompatible with the overall distribution are increased and decreased, respectively.
Furthermore, in complex biological data exhibiting a lot of noise, traditional statistical methods, such as the EM-algorithm, can struggle to cope with complex non-linear models when used to explore such data. In the Bayesian paradigm, on the other hand, all the unknown quantities are treated together in a consistent manner, to give fully probabilistic information on all unobserved variables, even their functions. Our method is based on a Chi-squared (χ 2)-mixture model for the gene-specific variances, with the number of components ranging from 1 to 5. In bioinformatics, as in many fields, mixture models have been fitted through the expectation-maximisation (EM) algorithm with different values of the number of mixture components [9, 17, 18]. In this paper, we present a Bayesian analysis of the variance mixture model, which we implement in the Bayesian software package, WinBUGS . The methodology is applied to a dataset on diffuse large B-cell lymphomas. The data contains expression levels of chronic lymphatic leukaemia (CLL) and diffuse large cell (DLCL) malignancies for the 9216 genes under study. The results of applying our model to the data-set are presented in the Results from the mixture model section. A discussion of the model and the results are in the Conclusions section.
It is hoped that the work presented in this paper will contribute to a large volume of current research work aimed at minimizing the risk of false positives in microarray experiments. The two extreme assumptions on the variance of gene expression data are presented in the Standard models for gene variance section and a description of the variance mixture model is contained in A mixture model for gene variance section.
We analyse the data described in , which investigates the classification of diffuse large B-cell lymphomas into distinct groups by gene expression profiling. Diffuse large B-cell lymphoma is an aggressive malignancy of mature B lymphocytes. It is estimated that, with an annual incidence of over 25,000 cases, it accounts for approximately 40% of all cases of non-Hodgkin's lymphoma. Currently, a combination of clinical parameters is used to provide an assessment of a patient's risk profile and to determine the most suitable clinical course of treatment. Whilst most patients initially respond well to chemotherapy, patients receiving the same diagnosis can have very different final outcomes in terms of remission achieved and their overall survival. It is suspected that the prognostic variables used are in fact proxies for the underlying cellular and molecular variation within diffuse large B-cell lymphomas. Their work considers whether gene expression profiling could subdivide this 'clinically heterogenous diagnostic category into molecularly distinct diseases that would possess more homogeneous clinical behaviours' . The microarrays used in this experiment were specially designed complementary DNA microarrays, called 'Lymphochips', which included those genes with a known or suspected role in processes that are important in immunology or cancer, and those genes known to be preferentially expressed in lymphoid cells. The profiling of gene expression included the three most prevalent adult lymphoid malignancies, in addition to profiling the gene expression in purified normal lymphocyte subpopulations under a range of activation conditions, in normal human tonsil and lymph node, and in a variety of lymphoma and leukaemia cell lines. From each experimental mRNA sample, a cDNA sample was prepared and labelled with red Cy5 dye. Furthermore, a corresponding reference cDNA sample, labelled with green Cy3 dye, was prepared from a pool of mRNAs isolated from nine different lymphoma cell lines. The labelled samples were combined and hybridized to the microarray.
Approximately 1.8-million measurements of gene expression were taken in all, across 96 normal and malignant lymphocyte samples, using 128 of the Lymphochip microarrays. To demonstrate the method proposed in this paper, we use a small subset, containing only eight slides from two conditions. Four slides quantify gene expression relating to chronic lymphatic leukaemia (CLL) malignancies and the other four to diffuse large cell lymphoma (DLCL) malignancies. Each slide contains measurements for 9216 genes. The red-to-green intensity ratio can be quantified for each gene and this reflects the relative abundance of mRNA in the experimental sample compared with the reference mRNA pool. By using a common reference sample, these fluorescent ratios can be considered a comparable measurement of the relative expression level of each gene across all of the samples. We want to compare different gene expression levels between the CLL and DLCL malignancy conditions.
where is gene i's log-ratio of observed intensity for condition c to the reference mRNA pool in replicate r c (i = 1, ..., 9216; c = 1 for CLL, 2 for DLCL; r 1 = 1, ..., 4; r 2 = 1, ..., 4); μic is the mean log-ratio; and is normally distributed with mean 0.
The mixture weights were estimated as part of the model, where they were assigned a Dirichlet prior distribution, (π 1, ..., π k ) ~ Dirichlet(1, ..., 1), where k is the number of mixture components. The χ 2 mixture scale parameters ψ j were assigned independent Gamma(0.01, 0.01) prior distributions. For each mixture model, three independent chains were run for 50,000 iterations. We discarded the first 20,000 iterations and used a combined sample of the remaining 60,000 iterations for posterior summaries.
Posterior mean (SD) for various variance mixture models
where m1 and m2 are sample sizes for condition 1 and 2, respectively.
Classification of genes under the Bayesian and EM mixture models, with four components
Top ten ranked genes by different variance models (Genes are listed by their codes)
Finally, other than setting , which is outright membership, we also considered a weighted variance estimate based on a four-component variance mixture model. We ranked genes using this weighted variance estimate to compare the ranking obtained from the other variance estimates (see Table 3). Eight genes coded 2143, 3181, 4069, 4323, 4532, 4586, 8076 and 8903 that appear in the top ten under the weighted variance model, also appear at least once in the top ten of the other variance models. These genes might be interesting, thus requiring further analysis and investigation.
We are aware that different assumptions for the prior distributions for both mixture weights and scales may give different results. We performed some limited sensitivity examination of the results to different specification of the priors. There were slight differences in the results, but the substantive conclusions were not affected.
We have presented a Bayesian variance mixture model for differential gene expression data. This model is a compromise between two extreme models: the too stringent constant gene variance and the overparameterised gene-specific variance models, which are both unrealistic assumptions. Our mixture variance model provides a more realistic and flexible estimate for the variance of gene expression data under limited replicates. We believe that in using the (weighted) latent class variances, estimated from a larger number of genes in each derived latent group, the p-values obtained are more accurate then either using a constant gene variance or gene-specific variance estimate.
In our example data, the results obtained from using our model are in close agreement to those obtained using EM algorithm implementation , which had been shown, using simulation studies, to be flexible and reliable in both true and false discovery rates. Our results are based on the assumptions of normally distributed log-ratios and a constant gene-variance between the two conditions. We are working on relaxing these conditions, in particular in using a long-tailed t-distribution as a robust alternative to allow for the possibility of gene-intensity measurement outliers. By varying the degrees of freedom, the t-distribution can also be used to investigate the sensitivity of the posterior results to changes in the prior for the gene intensity measurements.
We assume that the intensity level data is background corrected and normalised according to , using an arcsinh transformation based on a model for the dependence of the variance on the mean intensity levels with variance stabalizing properties. These are implemented in a variance stabilizing transformation vsn, the title function of the vsn package, part of the Bioconductor project for R . This results in generalised log-transformed expression intensity values. We develop the methodology for unpaired (indirect) data case, where two samples of interest are each co-hybridised with a reference sample. That is, each independent slide is a two-colour microarray experiment. The methodology is easily adapted to paired (direct comparison) data case.
whose usefulness depends on an accurate estimate of the standard error (SE) of i . Generally, there will be very few replications; thus the estimation of SE would be very unreliable and unstable. One solution is to discard genes with a small fold-change and very small standard deviations to avoid getting overoptimistic significant results.
There are two extreme cases to model this variance:
which has ν = m T + m C - 2 degrees of freedom. Under the null hypothesis, μ iT = μ iC , the standardised fold-change t i is distributed according to a t-distribution with ν degrees of freedom. This option has low power as it is based on very few replications. This makes estimation of the sample variances very unreliable and unstable, and this results in less powerful t-tests.
This has a large number of replicates, totalling (m T + m C )g. The degrees of freedom of the variance estimate are now ν = (m T + m C - 2)g, which makes the statistic t i behave like a standard Normal (0,1) variate. However, a constant variance over all the genes is too unrealistic an assumption, and it increases the risk of a false positive result for a gene with a larger variance. On the other hand, there is a high risk of missing out on a truly differentiated gene having a small variance, and a large differential effect.
where k is the number of latent components (classes); π = (π 1, ..., π k ) are the mixture proportions with π j being the probability that a gene belongs to latent class j (π 1 + ⋯ + π k = 1); ψ = (ψ 1, ..., ψ k ) are the component-specific χ 2 distribution scaling constants and g(w i |ψ j , ν) are χ 2 distributions with the scaled constant ψ j being specific to component j and ν = m T + m C - 2 is the common degrees of freedom, common for all components. The variance of the class j is = 1/ψ j , the variance of all the genes in class j.
In mixture modelling, it is convenient to formulate the model using a missing data problem, where each observation w i is assumed to arise from a specific but unknown component z i of the mixture. The model can be written in terms of the missing data z = z 1, ..., z g , otherwise known as allocation variables, which are assumed to be independent realizations of discrete random variables Z 1, ..., Z g with
Pr(Z i = j|π, ψ, ν) = π j (i = 1, ..., g; j = 1, ..., k).
The allocation of gene i to component j is based on the highest posterior probability π ij over all the components. Genes in the same latent class j share the same variance 1/ψ j , the mean variance of all the n j genes in component j, estimated from a larger number of replicates than purely those per gene (n 1 + … + n k = n). Parameters of the model can be estimated by the EM algorithm as shown in .
Quantities of interest, such as the posterior classification probabilities in (12), are estimated by plugging-in the point estimates and of π and ψ, respectively. Such plug-in estimates do not account for all the variability in estimating the model parameters and, as such, are more likely to underestimate the variance of the model parameters, which might inflate the significance levels. Aside from the problems associated with estimating the variability, the EM algorithm can sometimes have computational problems, not least in finding local maximum of the likelihood surface among several possible local maxima instead of the global maximum. In order to avoid the possibility of such problems, we propose a Bayesian hierarchical structure for the mixture model of the gene differential variance. We use an exchangeable gamma(ς, τ) prior on ψ i and a dirichlet(α 1, ..., α k ) prior on π. The hyperparameters ς and τ can be influential, and therefore in our full Bayesian analysis, these are not fixed, but given vague hyper-priors.
SOMM gratefully acknowledges the support of the Medical Research Council (MRC) (Grant Ref G0400908). REW was funded by the Engineering and Physical Sciences Research Council (EPSRC) on a postgraduate fellowship and MSG was funded by the Higher Education Funding Council of England (HEFCE). We thank two anonymous referees for their very helpful and constructive comments.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.