Skip to main content
  • Methodology article
  • Open access
  • Published:

Probabilistic estimation of microarray data reliability and underlying gene expression

Abstract

Background

The availability of high throughput methods for measurement of mRNA concentrations makes the reliability of conclusions drawn from the data and global quality control of samples and hybridization important issues. We address these issues by an information theoretic approach, applied to discretized expression values in replicated gene expression data.

Results

Our approach yields a quantitative measure of two important parameter classes: First, the probability P(σ|S) that a gene is in the biological state σ in a certain variety, given its observed expression S in the samples of that variety. Second, sample specific error probabilities which serve as consistency indicators of the measured samples of each variety. The method and its limitations are tested on gene expression data for developing murine B-cells and a t-test is used as reference. On a set of known genes it performs better than the t-test despite the crude discretization into only two expression levels. The consistency indicators, i.e. the error probabilities, correlate well with variations in the biological material and thus prove efficient.

Conclusions

The proposed method is effective in determining differential gene expression and sample reliability in replicated microarray data. Already at two discrete expression levels in each sample, it gives a good explanation of the data and is comparable to standard techniques.

Background

A broad variety of algorithms has been developed and used to extract biologically relevant information from gene expression data. Among others commonly used are visual inspection [1], hierarchical and k-means clustering [2], self organizing maps [3, 4] and singular value decomposition [5, 6]. These methods aim mainly at identifying predominant patterns and thus groups of "cooperating" genes based on the assumption that related genes have similar expression patterns.

Compared to the amount of work devoted to efficient methods to extract information from the data, somewhat less attention has been paid to the question of the reliability of the generated results. The ANOVA analysis [7] allows estimation, and thus elimination, of some systematic error sources. Bootstrapping cluster analysis estimates the stability of cluster assignments [8] based on artificial data-sets generated with ANOVA coefficients. Some authors also considered the question of how well a certain oligo [10] is suited to measure the mRNA expression level of the related gene.

Some work has gone towards the ambitious task of learning topological properties or qualitative features of the genetic regulatory network from expression profiles, see e.g. [11]. A major limiting factor in these attempts is the comparative sparseness of available data. It is therefore reasonable to consider reduced models, for example a Boolean representation of the gene activity. It is known that many biological properties, for instance stability and hysteresis, can be modeled by the dynamics of such reduced models [1214].

In this work we investigate the possibility of reducing complexity of gene expression data by discretizing the expression levels. The approach we present enables a new way of extracting biologically relevant information from the data in the following way: A biological variety, i.e. a biological system defined by the investigator, is represented by several samples which are subjected to gene expression analysis. If gene expression levels are discretized into n values, and the variety is represented by m samples, the number of observable expression states for a gene are limited to nm. These observed states S are modeled as being derived from a smaller number of underlying, biological states σ, through a measurement process. Rather than making static assignments S → σ we calculate conditional probabilities P(σ|S). The number of possible expression profiles for a gene over a set of varieties is limited and the probability of each expression profile is easily calculated. Since the model we use considers both the underlying biology and the measurement process it also generates a measure of sample coherence in each biological variety.

We demonstrate the feasibility of this approach for a binary discretization of gene expression. For the discretization step we use the absent/present classification provided by the Affymetrix software [9]. The outcome of our method on a data set covering gene expression in developing murine B-cells is compared to the results of a standard analysis. We show that even with the crude discretization into only two expression levels the method is competitive to statistical methods based on continuous expression levels.

Methods

The Model

A major step in the analysis of gene expression data is to separate the biological content of the data from measurement and sample specific errors. In other words given an observation, i.e. the expression values of a gene in several samples representing the same biological variety, one wants to conclude on the biological state σ, which generated the observation. This can be expressed as a conditional probability,

P(σ|S)     (1)

, that a gene is in a certain biological state σ given the corresponding observed state S. In the application on which we demonstrate the method we consider three different biological varieties: pro, pre, and mature B-cells. The samples in each variety are different cell lines arrested at the corresponding stage of development.

In this work we take an information theoretic point of view to estimate this probability: The information of interest, the state σ, is "transmitted" in a noisy measurement process and potentially distorted (Figure 1). Using Bayes' theorem, the desired conditional probability Eq. (1) can be expressed as:

Figure 1
figure 1

Schematic diagram illustrating the transition from underlying to observed distributions of states, in the case of m = 4 samples. The underlying distribution on the left hand side can be described by the probabilities for each underlying state, P1), P0), and Pr) (see text). This distribution is then distorted by sample specific errors, and , resulting in an experimentally observed distribution, depicted on the right hand side.

On the right hand side of this equation, P(S|σ) is the probability to observe state S if the underlying biological state is σ. In a sense, P(S|σ) describes the noise characteristic of the measurement process. In the following we will show how this conditional probability, and the other probabilities on the r.h.s. of Eq. (2) can be estimated.

Given a set of m samples representing the same biological variety, differences in the expression level of a gene between the samples can arise from two independent sources:

  1. 1.

    Random variation within the variety. This may be caused by temporal differences in response to the stimuli, slightly different environmental conditions, genotypic differences between samples, etc.

  2. 2.

    Sample specific errors. These are mainly caused by the measurement process, e.g. differences in the treatment of the mRNA, scratched arrays, and so on. However, outlier samples, cultured under considerably different conditions, also contribute to sample specific errors.

A separation of these two contributions is possible only with an appropriate model for the variation of gene expression between the samples. In the choice of model, one has considerable freedom within the bounds set by biological plausibility. A limiting factor on the biological model comes from the type and amount of available data. The data used in this work contains only four samples for each variety. For the model we propose this is the minimum number of samples required to estimate the model parameters.

In the discretization of gene expression levels, we use only two discrete values, 0 and 1, for the expression of a gene in a sample. This means that the number of observable states, S, in a variety consisting of m samples is 2m. With no measurement errors we could immediately conclude on the underlying biological state σ: the two cases, where all observations agree S = (1,...,1) and S = (0,...,0) can be mapped to the biological states σ1 and σ0 respectively, which describe "pure" states without variation. The remaining N - 2 observable states S, where the individual measurements disagree, correspond to biological states σ with random variation. For the application in our biological study with supposedly identical biological systems contributing to the observable states S, the exact pattern leading to contradicting observations does not carry any information, as long as we assume that there are no sample specific errors. Therefore, we subsummize all N - 2 possible observations as one biological state σr with a random variation.

The biological rationale for this model is given by the following example: If one considers a biological variety such as cells in the retina of the eye, then a certain number of crucial genes ought to be expressed in all samples. Such genes might include rhodopsin, a molecule that responds to light. In contrast, genes such as the hemoglobin family, which are typical of erythrocytes, ought not to be expressed in the retina. A third class of genes could be considered as independent of the system in the sense that their expression is not directly related to the biological system. Such genes may vary in expression both due to environmental and genetic differences between the samples.

The model discussed so far is depicted graphically in the left part of Figure 1, where a possible distribution of the relative frequencies of the three biological states is depicted, for the case of m - 4 samples. The distribution can be described by three numbers: the probabilities P1) and P0), which contribute to the frequencies of the states S = (1,1,1,1) and S = (0,0,0,0), and Pr) which contributes to both the frequency of mixed states and the two states above. Describing the mixed states with only one parameter Pr) implies that the biological variation is modeled evenly and identically distributed independently for each sample. In a second step, the measurement process with possible sample specific errors is modeled as statistically independent between samples. For each sample i, we define two parameters, and , denoted sample specific error probabilities.

To introduce the full formalism of our current model we start by considering a simple example, again for m - 4 samples. An observed state S, S ≡ (S1, S2, S3, S4) = (1,0,1,0), may be generated by the gene being in state σ1 with the probability:

or it may be generated by the gene being in state σ0 with the probability:

or it may be generated by the gene being in state σr with the probability:

With the briefer notation,

, where δ refers to the Kronecker delta (i.e. δj,k= 1 if j = k and 0 otherwise), we may express the distribution of observed states, in the general case of binary discretization with m samples, as:

Altogether the model uses 3 + 2 * m variables. These parameters P1), P0), Pr) and are estimated from the observed distribution of states (right side of Figure 1) by Levenberg-Marquardt [15] chi-square minimization of the unweighted error to the theoretical distribution Eq. (3). Using Eq. (2), and the parameters estimated as above, our belief that a gene belongs to the underlying states σ0, σ1, σr, given the 24 = 16 observable states S, can now be expressed as:

Once the probability that a gene is in a certain biological state Σi σ1, σ0, σr has been calculated for all varieties i = 1...v, one can calculate the probability that a gene exhibits a certain expression profile over a set of different varieties by taking the product

In this way, the probabilistic state analysis also generates a clustering: For a given expression profile over the varieties, e.g.

, we may extract those genes for which this expression profile is the most probable. In fact this is a "soft" clustering, in that an expression profile can belong to several clusters simultaneously with different probabilities. Moreover the genes clustered to a biologically interesting expression profile can be ranked by the probability of Eq. (4).

Experimental data preparation

All cells were grown in RPMI medium supplemented with 7.5% fetal calf serum, 10 mM HEPES, 2 mM pyruvate, 50 mM 2-mercaptoethanol and 50 mg gentamicin per ml (complete RPMI media) (all purchased from Life Technologies AB, T = E4by, Sweden) at 37 = B0C and 5% CO2. RNA was prepared using Trizol (GIBCO) and 7.5 = B5g of total RNA was annealed to a T7-oligo T primer by denaturation at 70 = B0C for 10 minutes followed by 10 minutes of incubation of the samples on ice. First strand synthesis was performed for 2 hours at 42 = B0C using 20 U of Superscript Reverse Transcriptase (GIBCO) in buffers and nucleotide mixes according to the manufacturers instructions. This was followed by a second strand synthesis for 2 hours at 16 = B0C, using RNAseH, E coli DNA polymerase I and E coli DNA ligase (all from GIBCO), according to the manufacturers instructions. The obtained double stranded cDNA was then blunted by the addition of 20 U of T4 DNA polymerase and incubation for 5 minutes at 16 = B0C. The material was then purified by Phenol:Cloroform:Isoamyl alcohol extraction followed by precipitation with NH4Ac and Ethanol. The cDNA was then used in an in vitro transcription reaction for 6 h at 37 = B0C using a T7 IVT kit and biotin labeled ribonucloetides. The obtained cRNA was purified from unincorporated nucleotides on a RNAeasy column (Qiagen). The eluted cRNA was then fragmented by incubation of the products for two hours in fragmentation buffer (40 mM Tris-acetate, pH 8.1, 100 mM KOAc, 150 mM MgOAc). 20 = B5g of the final fragmented cRNA was then hybridized to affymetrix chip U74Av2 (Affymetrix) in 200 = B5l hybridization buffer (100 mM MES-buffer, pH 6.6, 1 M NaCl, 20 mM EDTA, 0.01 Herring sperm DNA (100 = B5 g/ml) and Acetylated BSA (500 = B5 g/ml) in an Affymetrix Gene Chip Hybridization oven 320. The chip was then developed by the addition of FITC-streptavidin followed by washing using an Affymetrix Gene Chip Fluidics Station 400. Scanning was performed using a Hewlett Packard Gene Array Scanner.

Results

To evaluate the method we used both real and synthetic data. The experimental data was generated with Affymetrix microarrays for the study of differentiating murine B-cells at different stages in the differentiation process. In this publication the data is only used to demonstrate the feasibility of the proposed method. The biological implications of this study are published elsewhere [20].

Synthetic data and the effect of correlation

For synthetic data, generated with the model parameters P0) = 0.45, P1) = 0.35, Pr) = 0.2 and = = 0.02 for all samples i, parameter estimates are, as expected, given with low errors. The parameter values were chosen as typical values from the estimates on real data, see next section, and the result was verified for sample sizes m = 4, m = 5, and m = 6 (data not shown).

An assumption of simple model used to derive Eq. (3) is that randomly varying genes vary independently in the samples of a variety. Hence we investigated how severely this assumption influences the estimation of the model parameters.

To assess the influence of correlations between randomly varying genes we generated a data set consisting of four bits, i.e. samples, with the same parameters as above. In the random patterns a correlation was introduced between the third and fourth bit by changing the value of of the fourth to that of the third with a certain probability. We define this probability as the correlation factor. The correlation was introduced before distorting the patterns with error probabilities. We then plotted the mean error in the estimation of parameters over 500 runs of synthetically generated data for correlation factors in the range {0, 0.02,...,0.98}.

Figure 2 shows the error in the estimation of the parameters describing the underlying distribution. We notice that even for fully correlated patterns the estimation error is less than 20% of the correct values. The estimation of the probability for biologically varying genes is somewhat worse, for fully correlated patterns the error is almost 50%. For real data one can, however, expect a much smaller correlation. The average error in the estimates of the error probabilities, as seen in Figure 3, shows the expected behavior: The average error grows with the correlation for the uncorrelated samples, while the estimate for the correlated observations is almost unaffected. Intuitively, the model compensates for the correlation by increasing P1) and P0) as well as the error probabilities and lowering Pr). For correlation factors above 0.50, due to the compensation effect, the model deteriorates in explaining the data. This can be seen in the sum P0) + Pr) + P1) which initially drops from almost 1 to 0.99 as the correlation factor rises from 0 to 0.50 and then from 0.99 to 0.96 for correlation factors in the range 0.50 to 0.98 (data not shown). We hence conclude that it is reasonable not to impose the condition P0) + P1) + Pr) = 1 in the model, as this sum indicates if samples are strongly correlated in genes whose expression vary around the threshold of discretization.

Figure 2
figure 2

The average error in the estimation of the parameters P1), P0), Pr) are given as a function of correlation factor between the third and fourth bit. For correlation factors above 0.2 the error in Pr) rises considerably.

Figure 3
figure 3

The average error in the estimation of the error probabilities . For correlation factors above 0.2 and are notably raised. Patterns were these bits deviate from the other two are then not considered as random but rather caused by an error. This effect could only be avoided by introducing extra parameters for correlation between bits.

In summary, for not too large correlations in the biological variance the algorithm gives a good quantitative estimate of the model parameters. In the case of large correlations the qualitative picture given by the estimated parameters is still reliable.

Real data

Differentiating B-cells are characterized by phenotypic markers into different stages of development. Here we chose to study the expressional differences between three such stages; pro, pre and mature B-cells. For each of these three varieties we used four different cell lines arrested at the corresponding stage of development. Measurements we performed with Affymetrix array containing probesets for 12488 genes and ESTs on each sample. The discretization of expression levels was given by the Affymetrix GeneChip absent present calls [9].

Our algorithm was used to estimate the parameters P1), P0) and Pr), describing the biological distribution and the error probabilities (see Table 1). Theoretically, one expects the three biological probabilities to sum up to one. In our model, Eq. (3), we do not explicitly impose this condition. Nevertheless, the sum of the independently estimated parameters is close to one. This indicates that our model is a reasonable approximation of the biological system and the measurement process.

Table 1 Typical paramter values. Summary of the estimated parameter values for the B-cell data. P e refers to the set of error probabilities, i.e., .

The error probabilities from Eq. (3) can be used as a consistency index for the samples in a given variety. In the last variety (mature B-cells) the maximum error probability is notably higher. This effect is likely to be explained by the different anatomical origins of the cell lines representing this group. No such differences exist in the other groups since they all originate in the bone marrow which is the only anatomical site for B cell development in the adult animal [16]. In contrast, the mature B cell can reside in several other sites such as spleen, lymph-nodes and intestine which may affect the gene expression profile in these cells [17, 18]. With only four samples, it is not unlikely that these effects show up in the error probabilities and not only in the random variation parameter Pr).

Comparison to conventional t-test on known genes

To determine how well biologically relevant information can be extracted from the discretized data, we compare it with another statistical method based on continuous expression values. We use our method to identify differences in gene expression between two varieties in the following way. A gene that goes up between variety i and variety j is characterized by the states or or . Hence the belief that a gene goes up is given by the probability:

suppressing the conditional probabilities, P(·|S), for brevity. Similarly, the belief that a gene goes down is given by the probability:

Taking 1 - P(up) thus yields the Bayesian p-value of a gene going up. To answer the same question when working on continuous expression data one possibility is to employ a one sided two sample t-test in the Welch approximation of unknown variances in the varieties. This enables testing, for each gene, whether the mean of expression is higher or lower in one variety than in another. For comparison of these two approaches we selected a set of genes based on their well documented expression pattern and biological functions in the developing B lymphocyte [16, 19]. Several of these are functionally linked since they participate directly in somatic DNA rearrangement events occurring specifically at the pre-B cell stage or participate in the regulation of genes involved in this process and thus display restricted expression patterns (pre-B specific). A second set of genes were selected based on their expression in cells that are either committed to the B lineage (B-lineage specific genes, in pre-B and B-cells) or non committed to this developmental pathway (Not in B-lineage, expressed in pro-B cells) [21].

The result of this comparison is presented in Table 2. For 14 out of the 22 genes the two methods completely agree. Out of these 14 only one (Mb1) does not match the expected target profile. For the other 8 genes, where the two methods yield different results, the probabilistic state analysis gives the expected answer in 5 cases, which should be compared to the two cases, where the t-test gives the right answer. In one case (rag-1), neither of the two methods gives the expected result.

Table 2 t-test vs. probabilistic analysis of gene expression levels. The three groups I, II and III indicate the expressional changes between Pro-B to Pre-B, Pre-B to Mature-B, and Pro-B to Mature-B respectively. U stands for accepting the hypothesis up, D for down, and S (stable) if no hypothesis could be accepted on the 95% confidence level.

For the subset of genes considered here, our method has an advantage of 5 : 2 in giving the correct (i.e. expected) expression pattern. However, the number of samples is not big enough to draw firm conclusions from this result.

Conclusions

The method we have presented serves several purposes:

  1. 1.

    It gives a measure of the biological variation of the genes' expression in different varieties.

  2. 2.

    It estimates each hybridizations' global error probabilities. These parameters are very useful as they serve as quality/consistency indicators of the samples of each variety.

  3. 3.

    Given the parameters above, it estimates the probability of a gene belonging to each of the three groups σ0, σr and σ1. These probabilities in turn indicate weather the gene is likely to be below, fluctuating around or above the threshold of discretization.

  4. 4.

    Clustering of genes to expression profiles over a set of different varieties is achieved with Eq. (4). The probability, i.e. belief, that a gene belongs to a certain cluster is exactly quantified.

This novel approach is proven valuable for quantifying both data reliability and underlying gene expression in microarray experiments. Our method has been successfully applied in two different projects [22, 20].

References

  1. Cho RJ, Campbell MJ, Winzeler EA, Steinmetz L, Conway A, Wodicka L, Wolfsberg TG, Gabrielian AE, Landsman D, Lockhart DJ: Genome-Wide transcriptional analysis of the mitotic cell cycle. Mol Cel 1998, 2: 65–73.

    Article  CAS  Google Scholar 

  2. Eisen MB, Spellman PT, Brown PO, Botstein D: Cluster Analysis of Genome Wide Expression Patterns. PNAS 1998, 95: 14863–68. 10.1073/pnas.95.25.14863

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  3. Tamayo P, Slonim D, Mesirov J, Zhu Q, Kitareewan S, Dimitrovsky E, Lander E, Golub TR: Interpreting Gene Expression with Self-Organizing Maps. PNAS 1999, 96: 2907–2912. 10.1073/pnas.96.6.2907

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  4. Törönen P, Kolehmainen M, Wong G, Castren E: Analysis and Visualization Of Gene Expression Data Using Self-Organizing Maps. FEBS Lett 1999, 451: 142–146. 10.1016/S0014-5793(99)00524-4

    Article  PubMed  Google Scholar 

  5. Holter NS, Mitra M, Maritan A, Cieplak M, Banavar JR, Fedoroff N: Fundamental patterns underlying gene expression profiles: Simplicity from complexity. PNAS 2000, 97: 8409–8414. 10.1073/pnas.150242097

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  6. Alter O, Brown P, Botstein D: Singular value decomposition for genome-wide expression data processing and modeling. PNAS 2000, 97: 10101–10106. 10.1073/pnas.97.18.10101

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  7. Kerr MK, Martin M, Churchill GA: Analysis of Variance for Gene Expression Microarray Data. J Comp Biol 2000, 7: 819–837. 10.1089/10665270050514954

    Article  CAS  Google Scholar 

  8. Kerr MK, Churchill GA: Bootstrapping cluster analysis: Assessing the reliability of conclusions from micro-array experiments. PNAS 2001, 98: 9061–8965.

    Google Scholar 

  9. Affymetrix technical notes[http://www.affymetrix.com]

  10. Li F, Stormo GD: Selection of optimal DNA oligos for gene expression arrays. Bioinformatics 2001, 17: 1067–1076. 10.1093/bioinformatics/17.11.1067

    Article  CAS  PubMed  Google Scholar 

  11. Friedman N, Linial M, Nachman I, Peér D: Using Bayesian Networks to Analyze Expression Data. J Comp Biol 2000, 7: 601–620. 10.1089/106652700750050961

    Article  CAS  Google Scholar 

  12. Kauffman SA: Metabolic Stability and Epigenesis in Randomly Constructed Genetic Nets. J Theor Biol 1969, 22: 437.

    Article  CAS  PubMed  Google Scholar 

  13. Kauffman SA: The Origins of Order. Oxford University Press 1993.

    Google Scholar 

  14. Huang S: Gene expression profiling, genetic networks and cellular states: an integrating concept of tumor genesis and drug discovery. J Mol Med 1999, 77: 469–480. 10.1007/s001099900023

    Article  CAS  PubMed  Google Scholar 

  15. Marquardt DW: Journal of the Society for Industrial and applied Mathematics 11: 431–441.

  16. Gia P, ten Boekel E, Rolink AG, Melchers F: B-cell development: a comparison between mouse and man. Immunol Today 1998, 19: 480–5. 10.1016/S0167-5699(98)01330-9

    Article  Google Scholar 

  17. Rolink AG, Melchers F, Andersson J: The transition from immature to mature B cells. Curr Top Microbiol Immunol 1999, 246: 39–43.

    CAS  PubMed  Google Scholar 

  18. Rolink AG, ten Boekel E, Yamagami T, Ceredig R, Andersson J, Melchers F: B cell development in the mouse from early progenitors to mature B cells. Immunol Lett 1999, 68: 89–93. 10.1016/S0165-2478(99)00035-8

    Article  CAS  PubMed  Google Scholar 

  19. Liberg D, Sigvardsson M: Transcriptional regulation in B cell differentiation. Crit Rev Immunol 1999, 19: 127–53.

    CAS  PubMed  Google Scholar 

  20. Tsapogas P, Breslin T, Bilke S, Lagergren A, Mnsson R, Liberg D, Peterson C, Sigvardsson M: RNA analysis of B-cell lines arrested at defined stages of differentiation allows for an approximation of gene expression patterns during B cell development. J Leukoc Biol, in press.

  21. Rolink AG, Schaniel C, Busslinger M, Nutt SL, Melchers F: Fidelity and infidelity in commitment to B-lymphocyte lineage development. Immunol Rev 2000, 175: 104–11. 10.1034/j.1600-065X.2000.017512.x

    Article  CAS  PubMed  Google Scholar 

  22. Högerkorp CM, Bilke S, Breslin T, Ingvarsson S, Borrebaeck C: CD44-stimulated human B cells express transcripts specifically involved in immunomodulation and inflammation as analyzed by DNA microarrays. Blood 2003, 101: 2307–13. 10.1182/blood-2002-06-1837

    Article  PubMed  Google Scholar 

Download references

Acknowledgments

SB and TB are supported by Knut and Alice Wallenberg Foundation through the SWEGENE consortium. MS is supported by the Swedish science council and the A. Österlund foundation. The authors also extend their gratitude to Carsten Peterson, Patrik Edén, Jari Häkkinen and Markus Rignér, for fruitful discussions and careful proofreading of the manuscript.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Sven Bilke.

Additional information

Authors' Contributions

SB and TB have contributed to the development of the model and implemented the algorithms. MS has contributed the experimental data and biological expertise.

Authors’ original submitted files for images

Rights and permissions

Reprints and permissions

About this article

Cite this article

Bilke, S., Breslin, T. & Sigvardsson, M. Probabilistic estimation of microarray data reliability and underlying gene expression. BMC Bioinformatics 4, 40 (2003). https://doi.org/10.1186/1471-2105-4-40

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2105-4-40

Keywords