<title>Introduction</title> <p>In examining genes, either individually or in system-wide characterizations, it is useful to generalize its “state”. For example, a gene's <italic>Present/Absent</italic> call is a common dimension of the reported results of gene-expression microarray experiments. Such calls tag each probe set in the microarray with a determination of whether the probe set is expressed (Present) or unexpressed (Absent) in the sampled tissue ##REF##12490443##[1]##., <italic>Present/Absent</italic> calls are often used in filtering out false positives from the large collection of probes on an expression array. The most commonly used approach to making such calls is the MAS5 algorithm ##REF##12490443##[1]##, part of the Affymetrix™ collection of software tools ##REF##15846360##[2]##. While some recent experimental findings support the use of the MAS5 algorithm ##REF##15693945##[3]##, MAS5 has some significant shortcomings. First, MAS5 does not provide the user with a statistical gauge of the basic claim behind the <italic>Present/Absent</italic> call. Second, MAS5 does not compare calls across multiple samples. Finally, because MAS5 does not operate on adjusted readings, it cannot benefit from the increasingly sophisticated techniques for adjusting gene expression readings (e.g. RMA ##REF##12582260##[4]## and others ##UREF##0##[5]##; see ##REF##16113766##[6]## for a comparison of techniques)</p> <p>Conceptually it is understood that the classification of genes into alternative states is a simplification of much greater complexity patterns of gene behaviour and action. However, empiric evaluation of the observed data finds that gene expression patterns commonly can fit one of two alternative expression level distributions. Moreover, such simplification has proven valuable in other research domains. For example the simplification that abstracts digital logic from the underlying continuous flow of electrons in integrated circuits has enabled the design of devices of staggeringly complex functionality ##REF##12353013##[7]##.</p> <p>We describe here a method that makes use of quantitative expression level readings. It is important to stress that the method is not a pre-processing step, like background adjustment for noise, but rather a post-processing step that makes use of the noise-adjusted readings. In the specific examples presented here, we make use of RMA-adjusted expression levels ##REF##12582260##[4]## from Affymetrix microarrays, but the input could be raw or adjusted values from any platform. Using the expression levels, we build a statistical model of expression for each probe set, based on an assumed bimodal distribution, that accounts for the two states of an expressed gene: Up and Down.</p> <p>The inputs to the statistical model are the probe-set expression levels from multi-sample experiments. For a specific probe set, we gather expression levels from the cohort of samples for. For example, in a set of experiment involving 100 patient samples and 100 control samples, we obtain 200 data points for the single probe set (see ##FIG##0##Fig. 1##). We then use the data points from the single probe set to infer two gamma distributions, one distribution representing the <italic>Down</italic> state and one representing the <italic>Up</italic> state. Our choice of gamma distributions comes from the distribution flexibility in containing the two distribution shapes we required. Such mixture models have been successfully applied to other problems in biology (e,g, ##UREF##1##[8]##–##REF##16204095##[10]##).</p> <p>A gamma distribution has the general form:where <italic>a</italic> is often called the shape parameter and <italic>b</italic> the scaling parameter. The general form of a gamma distribution is beyond the scope of this paper (see ##UREF##2##[11]## as a reference). For smaller values of <italic>a</italic>, the gamma distribution takes an exponential-like form, with a continuous decay that starts at zero; for larger values of <italic>a</italic>, the distribution takes a form similar to the normal distribution, with a mean of <italic>ab</italic> and variance of <italic>a</italic>−<italic>b</italic>
<sup>2</sup>.</p> <p>By combining single probe data across multiple samples, we consider the entire population of probe expression values (gene values) as derived from a single distribution. That single distribution is in fact the mixture of two gamma distributions – one distribution for the Up state and one distribution for the Down state. We represent the resulting model with six parameters: <italic>a<sub>u</sub></italic>, the shape parameter for the Up distribution; <italic>b<sub>u</sub></italic>, the scale parameter for the Up distribution; <italic>a<sub>d</sub></italic>, the shape parameter for the Down distribution; <italic>b<sub>d</sub></italic>, the scale parameter for the Down distribution, and <italic>η<sub>u</sub></italic>, <italic>η<sub>d</sub></italic>, the mixture coefficients that give the relation between the two distributions in the final mixture. We determine values for the different parameters using an Expectation-Maximization (EM) algorithm (see <xref ref-type="sec" rid="s4">Methods</xref>), the output of which are the six defining parameters. Upon completion of processing using the gamma mixture (GM) algorithm, we are able to calculate, given a specific expression value, the probability this expression value represents a gene in the <italic>Up</italic> (or <italic>Down</italic>) state.</p>
<title>Methods</title> <title>EM algorithm</title> <p>For each probed gene, the algorithm determine six parameters that define, together, the coefficients for each of the distribution (Up, Down) and the mixture coefficients between the two distributions. We call the set of different parameters <italic>θ</italic>
Where <italic>a<sub>u</sub></italic>, <italic>b<sub>u</sub></italic> determine the coefficient of the Gamma distribution that describes the Up gene state; <italic>a<sub>d</sub></italic>, <italic>b<sub>d</sub></italic> determine the coefficients of the Down state and <italic>η<sub>u</sub></italic>, <italic>η<sub>d</sub></italic> determine the mixture coefficient (and <italic>η<sub>u</sub></italic>, +<italic>η<sub>d</sub></italic> = 1).</p> <p>The algorithm is iterates over the different function, so that every iteration improves the estimate of the coefficients. In ##UREF##4##[14]## you can see the general proof of the EM algorithm, according to which, it is sufficient to find maximas for the function <italic>Q</italic>, defined as:Where <italic>θ</italic> is defined previously as the collection of parameters. <italic>θ</italic>
<sup>0</sup> stands for the set of parameters at a previous iteration and the index <italic>i</italic> goes over the two different function in the mixture and the index <italic>t</italic> goes over available data points.Finding maximas of <italic>Q</italic> replaces (the harder task of) finding maximas for the original function. To find maximas for <italic>Q</italic>, we differentiate it with respect to the model parameters and compare to zero. First according to <italic>b<sub>i</sub></italic>:
And then according to <italic>a<sub>i</sub></italic>:Where Ψ(<italic>x</italic>) is the psi function .</p> <p>Using a Lagrange multiplier to incorporate the constrainwe have to maximize the target functionwith respect to the <italic>η<sub>i</sub></italic>, we deriveand obtainWe solve this numerically (using Matlab®) in every iterative step, until we reach some predefined convergence criterion</p> <title>Gene expression spike-in data</title> <p>We used Affymetrix's deposited gene expression data for a SpikeIn experiment, as it is available from ##UREF##0##[5]##.</p> <p>Adjustments of the data were made using the RMAExpress tool ##UREF##5##[15]## over original CEL files. Affymetrix's Present/Absent call (MAS5 calls) were made using Affymetrix's GCOS tools ##REF##12490443##[1]##.</p> <p>Other sources of gene expression data:</p> <p>Data for the set of Bittner et. al were made available from data made publicly available by the Expression Project for Oncology, an International Genomics Consortium public/private initiative ##UREF##3##[13]##. Data from Miller el. al ##REF##16141321##[12]## has been obtained from the Gene Expression Omnibus ##REF##11752295##[16]##.</p>
<title>Results</title> <p>To compare the consistency of MAS5 calls with the consistency of GM calls, we used the publicly available results of a spike-in experiment ##UREF##0##[5]##. In this experiment, the researchers assayed samples that were identical except for controlled differences in the RNA of 42 transcripts. Except for the 42 transcripts whose levels were systematically manipulated, each probe would be expected to have the same <italic>Present/Absent</italic> call across the experiment and to have the same <italic>Up/Down</italic> call across the experiment. We measure the success of the two algorithms by their consistency over the cohort of sample. A perfect score for an algorithm would mean that the algorithm succeeded in finding identical <italic>Present/Absent</italic> or <italic>Up/Down</italic> call for each of the genes across the experiment.</p> <p>Of the 22,283 probes examined in the experiment, the MAS5 algorithm was consistent in assigning the same <italic>Present/Absent</italic> call, across all samples, for each of 17,004 probes; the remaining 5278 probes were assigned inconsistent calls by MAS5. In contrast, the GM algorithm consistently assigned the same <italic>Up/Down</italic> call for each of 19,923 probes and gave inconsistent calls for the remaining 2359 probes. Thus the GM algorithm showed an improvement of 55% in consistency.</p> <p>To examine the performance of the algorithms on data with natural biological variation, we turned to other publicly available studies. One such study, Miller et. al. ##REF##16141321##[12]##, provides U133-A/B data on 251 primary invasive breast tumor samples. We are especially interested in the ability of the MAS5 and GM algorithms to make calls that are consistent with (RMA-adjusted) expression levels. That is, we expect an <italic>Absent</italic> or <italic>Down</italic> call to correlate with low levels of expression and a <italic>Present</italic> or <italic>Up</italic> call to correlate with high readings. ##FIG##1##Figure 2## shows an example, probe set ‘206378_at’ (which represents the gene SCGB2A2), where these expectations are confounded. Panel (a) shows a simple histogram of expression levels from the probe set, across all samples; Panel (b) shows the derived probability distribution, based on the Gamma Mixture hypothesis; and Panel (c) plots the probability of being in an <italic>Up</italic> state, as a function of the expression level. As the figure shows, the <italic>Up/Down</italic> classifications produced by GM algorithm correlate well with expression values, across the range of expression values. The MAS5 algorithm, on the other hand, toggles between <italic>Present/Absent</italic> calls quite sporadically in the expression range. To compare the MAS5 calls and the GM calculated probability over a large set of samples, we made use of data from ##UREF##3##[13]##, following the procedure described in <xref ref-type="sec" rid="s4">Methods</xref> to obtain <italic>Present/Absent</italic> and <italic>Up/Down</italic> readings. ##FIG##2##Figure 3## shows the different readings. Panel (a), (b) and (c), as before, show the expression distributions of the probe set. In panel (d) a zoom-in view of the transition area of panel(c), shows the details of decision of the shift between the <italic>Down</italic> and <italic>Up</italic> status. In Panel (e), we can see the differences between decision based on the GM algorithm and the MAS5 algorithm, where low levels of expression values are toggled between present and absent calls made by MAS5 and, on the other hand, have low probability to be in the Up state (or high probability to be in the Down state).</p>
<title>Discussion</title> <p>We have described a new approach to determining the expression state of a gene. Like the Affymetrix MAS5 algorithm, our method is a two-state classifier. In contrast to the MAS5 algorithm, our method takes account of the underlying distribution of expression values in a set of samples. In particular, our method assumes a two-state distribution of gene expression that can be captured by mixed gamma distributions. As we have shown, our technique yields more stable calls than MAS5. In a set of biological replicates, MAS5 produced inconsistent calls for twice as many probe sets as the GM algorithm. Further, in a set of samples showing normal biological variation, the GM algorithm yielded calls that had better correlations with RMA-adjusted expression levels than the MAS5 calls.</p> <p>
##TAB##0##Table 1## shows the strengths and weaknesses of the two approaches. One of the main differences – both a weakness and a strength – is the fact that the MAS5 algorithm is applied to individual samples, while the GM algorithm is applied to a set of samples. On the one hand, this allows MAS5 to be applied to individual samples from arbitrary experiments; on the other hand, MAS5 cannot take advantage of the statistical power in the multi-sample joined population. Second, the GM approach is not limited to readings from Affymetrix platforms, but may be applied to any values that represent gene expression or protein abundance. Third, since MAS5 is applied to raw values, it cannot take advantage of new adjustment techniques such as RMA; in contrast, the GM algorithm can be applied to raw or adjusted values. . Fourth, while MAS5 does give a p-value for the Present/Absent decision, this p-value cannot be interpreted as a distance from population and does not convey biological information.</p> <p>In summary, we believe our approach to be a general and powerful way to fit gene expression data to a two-state model. We consider the GM call to be a true, scale free, normalization that is entirely platform-independent, applicable to any gene expression. While applied to gene expression microarrays that measure RNA abundance, this method is applicable to any quantitative measure of individual gene state.</p>
<p>Conceived and designed the experiments: SE LC KHB. Performed the experiments: SE. Analyzed the data: SE LC CFS. Wrote the paper: SE CFS KHB.</p> <p>We introduce a novel technique to determine the expression state of a gene from quantitative information measuring its expression. Adopting a productive abstraction from current thinking in molecular biology, we consider two expression states for a gene - <italic>Up</italic> or <italic>Down</italic>. We determine this state by using a statistical model that assumes the data behaves as a combination of two biological distributions. Given a cohort of hybridizations, our algorithm predicts, for the single reading, the probability of each gene's being in an <italic>Up</italic> or a <italic>Down</italic> state in each hybridization. Using a series of publicly available gene expression data sets, we demonstrate that our algorithm outperforms the prevalent algorithm. We also show that our algorithm can be used in conjunction with expression adjustment techniques to produce a more biologically sound gene-state call. The technique we present here enables a routine update, where the continuously evolving expression level adjustments feed into gene-state calculations. The technique can be applied in almost any multi-sample gene expression experiment, and holds equal promise for protein abundance experiments.</p>
<fig id="pone-0002901-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002901.g001</object-id><label>Figure 1</label><caption><title>Data handling in the mixture algorithm.</title><p>First, gene expression data from a set of gene expression experiments is collected. The matrix in the figure shows rows probe set, where every row is a single probe set, and every column is a different hybridization experiment. This could be, for example, Affymetrix microarray experiments, where each column is a different patient. We then look at the data probe-by-probe. For example, we follow probe “a” in the figure and look at the expression levels for this probe, across all samples in the set of gene expression experiments. Each probe will have data from the entire collection of experiments. For the specific probe “a”, we fit the set of expression measurements into two gamma distributions, one representing the “down” state and one representing the “up” state. Each data point is then computationally associated with a probability of being either under the first Gamma distribution (which would mean the gene associated with the probe, for the specific sample, is at a “down” state) or with the second Gamma distribution (which would mean the gene associated with the probe, for the specific sample, is at an “up” state). We iterate the procedure across the entire probe-set, to tag every gene across the microarray with its probability of being “up” or “down”.</p></caption></fig> <fig id="pone-0002901-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002901.g002</object-id><label>Figure 2</label><caption><title>An in-depth look into data from a single probe over a collection of 251 hybridizations of breast cancer samples.</title><p>(a) Displays a histogram of gene expression data for a single probe across the collection of 251 samples. The x-axis corresponds with levels of expression and the y-axis is a count of expression level for the specific bin. As the panel shows, many of the hybridizations show a level of expression close to zero. This is visible through a large collection (large count) of gene expression measurements close to zero expression levels. On the other hand, many of the probes show expression levels that spread across the entire span of expression levels from zero to 16,000. This is well-fitted into a gamma mixture model that assumed a dual behavior for the gene. Panel (b) shows a plot of the approximated distribution functions across the entire gene expression range for the surveyed probe. The two lines plotted have been calculated by fitting gene expression into the assumption of a binary state distribution, with each distribution modeled by a Gamma-like behavior. Panel (c) gives the probability of being in an “up” state as a function of gene expression for the specific probe surveyed. As panels (b) and (c) overlap, we demonstrate how changes in expression levels in (c) associate that level with the curves of panel (b). The higher an expression level in (c) is, the more probable it is to be affiliated with the red curve (“up” curve) of panel (b). The lower the expression level is in (c), it is more probable to affiliate it with the blue curve (“down” curve) of panel (b), and with a “down” state. Being in a “down” state is the reciprocal of being in an “up” state, which gives a probability of zero for being “up”.</p></caption></fig> <fig id="pone-0002901-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002901.g003</object-id><label>Figure 3</label><caption><title>Similar to matching panels (a) (b) and (c) in ##FIG##1##Figure 2##, panels (a), (b) and (c) of this figure show the histogram of gene expression, two resulting gamma curves, and the probability of being in an “up” state for a specific single probe set out of the collection of probes set in a collection of sample.</title><p>Panel (d) is a zoom into the highlighted part of panel (c). Panel (e) is a zoom into the highlighted part of panel (d). Panel (d) shows the gradual probabilistic transition from being associated with a “down” state to being in an “up” state. The transition correlates well with gene expression and demonstrates the sensitivity of the approach to changes in gene expression. In panel (d) and in zoom-in panel (e), we also highlight the Present/Absent calls made by MAS5. Especially in the panel (e), it is easy to see how MAS5 Present/Absent calls toggle with growing levels of expression, despite an expected plateau, the MAS5 algorithm stabilizes on a Present call at much higher levels of expression and makes an Absent call for gene expression level as high as 500, while giving a Present call to expression level of 50. Use the toggled calls puts the user in danger of associating very different of expression level with very different states of a gene.</p></caption></fig>
<table-wrap id="pone-0002901-t001" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002901.t001</object-id><label>Table 1</label><alternatives><table frame="hsides" rules="groups"><colgroup span="1"><col align="left" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/></colgroup><thead><tr><td align="left" rowspan="1" colspan="1">Technique</td><td align="left" rowspan="1" colspan="1">Potential use</td><td align="left" rowspan="1" colspan="1">Ability to incorporate new adjustment techniques</td><td align="left" rowspan="1" colspan="1">Statistical significance to results</td><td align="left" rowspan="1" colspan="1">Platform dependent</td></tr></thead><tbody><tr><td align="left" rowspan="1" colspan="1">MAS5</td><td align="left" rowspan="1" colspan="1">Affymetrix-based. Can be used on a single samples.</td><td align="left" rowspan="1" colspan="1">No</td><td align="left" rowspan="1" colspan="1">p-value</td><td align="left" rowspan="1" colspan="1">Yes</td></tr><tr><td align="left" rowspan="1" colspan="1">GM (gamma mixture)</td><td align="left" rowspan="1" colspan="1">Any gene expression reading. Can be used only in multi-sample experiments</td><td align="left" rowspan="1" colspan="1">Yes.</td><td align="left" rowspan="1" colspan="1">Yes</td><td align="left" rowspan="1" colspan="1">No</td></tr></tbody></table></alternatives></table-wrap>
<disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <inline-formula></inline-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula>
<fn-group><fn fn-type="COI-statement"><p><bold>Competing Interests: </bold>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p><bold>Funding: </bold>The authors have no support or funding to report.</p></fn></fn-group>
<graphic xlink:href="pone.0002901.g001"/> <graphic xlink:href="pone.0002901.e001.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.g002"/> <graphic xlink:href="pone.0002901.g003"/> <graphic id="pone-0002901-t001-1" xlink:href="pone.0002901.t001"/> <graphic xlink:href="pone.0002901.e002.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e003.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e004.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e005.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e006.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e007.jpg" mimetype="image" position="float"/> <inline-graphic xlink:href="pone.0002901.e008.jpg" mimetype="image"/> <graphic xlink:href="pone.0002901.e009.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e010.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e011.jpg" mimetype="image" position="float"/> <graphic xlink:href="pone.0002901.e012.jpg" mimetype="image" position="float"/>
[{"label": ["5"], "element-citation": ["\n"], "surname": ["Irizarry", "Jaffee"], "given-names": ["RA", "HA"], "year": ["2007"], "article-title": ["Affycomp II"]}, {"label": ["8"], "element-citation": ["\n"], "surname": ["Lee", "Wang", "Yau", "McLachlan", "Ng"], "given-names": ["AH", "K", "KK", "GJ", "SK"], "year": ["2007"], "article-title": ["Maternity Length of Stay Modelling by Gamma Mixture Regression with Random Effects."], "source": ["Biom J"]}, {"label": ["11"], "element-citation": ["\n"], "surname": ["Hogg", "McKean", "Craig"], "given-names": ["RV", "JW", "AT"], "year": ["2005"], "source": ["Introduction to mathematical statistics"], "publisher-loc": ["Upper Saddle River, N.J."], "publisher-name": ["Pearson Education"], "fpage": ["xiii, 704"]}, {"label": ["13"], "element-citation": ["\n"], "surname": ["Bittner"], "given-names": ["M"], "year": ["2006"], "article-title": ["International Genomics Consortium (IGC)"]}, {"label": ["14"], "element-citation": ["\n"], "surname": ["Duda", "Hart", "Stork"], "given-names": ["RO", "PE", "DG"], "year": ["2001"], "source": ["Pattern classification"], "publisher-loc": ["New York"], "publisher-name": ["Wiley"], "fpage": ["xx, 654"]}, {"label": ["15"], "element-citation": ["\n"], "surname": ["Bolstad"], "given-names": ["B"], "year": ["2007"], "article-title": ["RMAExpress. 1.0 ed"]}]
PLoS One. 2008 Aug 6; 3(8):e2901