A New Distribution Family for Microarray Data

Abstract

The traditional approach with microarray data has been to apply transformations that approximately normalize them, with the drawback of losing the original scale. The alternative standpoint taken here is to search for models that fit the data, characterized by the presence of negative values, preserving their scale; one advantage of this strategy is that it facilitates a direct interpretation of the results. A new family of distributions named gpower-normal indexed by $p\in \mathbb{R}$ is introduced and it is proven that these variables become normal or truncated normal when a suitable gpower transformation is applied. Expressions are given for moments and quantiles, in terms of the truncated normal density. This new family can be used to model asymmetric data that include non-positive values, as required for microarray analysis. Moreover, it has been proven that the gpower-normal family is a special case of pseudo-dispersion models, inheriting all the good properties of these models, such as asymptotic normality for small variances. A combined maximum likelihood method is proposed to estimate the model parameters, and it is applied to microarray and contamination data. R codes are available from the authors upon request.

1. Introduction

While analysing microarray intensity measurements, it is usual to find asymmetric distributions with some negative values and the purpose of this article is to model data with these characteristics.

The traditional approach with microarray data has been to apply transformations that approximately normalize them, with the drawback of losing the original scale. The initial transformation applied was $lo{g}_2$; it allows working with log-ratios which have a simple and intuitive meaning for biologists (see for example [1,2]). This transformation usually works well for high values but not for zero, and low ones. Further, it cannot be applied to negative values. To avoid these drawbacks, [3,4] suggested the generalized logarithm transformation (glog), that allows negative values and this transformation is obtained from a multiplicative–additive linear error model for the data, through a Taylor approximation.

On the other hand, the glog transformation usually works well for low values but it is too severe for high ones. The next improvement was introduced by [5], who defined transformations on a real supported data family named generalized power transformations (gpower):

$$gpower\left(Y;p\right)=\left\{\begin{array}{cc}\frac{{\left(Y+\sqrt{Y^2+1}\right)}^p-1}{p}\hfill & \mathrm{if}p\ne 0\hfill \\ ln\left(Y+\sqrt{Y^2+1}\right)\hfill & \mathrm{if}p=0\hfill \end{array}\right.$$

(1)

The gpower transformations extend the glog transformation continuously, in the same sense as the Box-Cox family [6] extends the natural logarithm.

A related problem presented by microarray data when they are log transformed is the intensity-dependent biases observed in the MA (Minus Average) plots (see [7]). These plots display pairwise comparisons of log intensities ($Int1$ and $Int2$) between microarrays. More specifically the vertical axis gives the log ratio or difference of the intensities in a log scale $M=lo{g}_2(Int1/Int2)=lo{g}_2\left(Int1\right)-lo{g}_2\left(Int2\right)$ and the horizontal axis is the average log intensities $A=1/2lo{g}_2(Int1\times Int2)=1/2(lo{g}_2\left(Int1\right)+lo{g}_2\left(Int2\right))$. Even when intensities from a controlled spike-in experiment are expected to show a horizontal cloud centred on the $M=0$ axis from low to high intensities, non-horizontal structures are apparent. Several approaches have been taken to account for the observed biases (see for example [2,8,9,10]) , explain and correct the observed MA plot intensity-dependent biases through linear transformations.

To avoid transformations, the alternative standpoint taken here is to search for models that fit the data, preserving their scale. One advantage of this strategy is that it facilitates a direct interpretation of the results. In this direction, [11] showed that data that become normal after a glog transformation belong to what they called the glog-normal distribution family.

In this paper, we extend their results by characterizing those distributions that become normal (or truncated normal) after a gpower transformation. We introduce the gpower-normal family; this family of distributions should be fitted to gene intensities that have been previously calibrated with an affine transformation, according to the Bengtsson and Hössjer proposal [10].

For positive data, a study has been carried out by [12], who analysed the power normal family, related to Box-Cox transformation. The improvement of the gpower-normal family over the glog-normal family is that it can account for lighter asymmetries where the glog transformation is too strong (see [1,5]); also, it has support on the whole real line, not being restricted to positive values.

In Section 2, we describe the development of a new probability model to be used as a statistical tool for microarrays data analysis. Gpower-normal models are defined and their main properties are demonstrated; their relation with pseudo-dispersion models is studied and expressions for the moments and quantiles are obtained. Then, a combined maximum likelihood method is described to obtain estimators of the parameters and it is applied in subsequent sections. For the purpose of illustration, in Section 3, we show several density functions, their quantiles and real data applications. Finally, discussion and conclusions are presented in Section 4 and Section 5 respectively. We have placed some proofs in the Appendices so as not to break the flow of the narrative.

2. Methods

The methodological idea within this section is to obtain a model for microarray data in their original scale. With this purpose, we describe a new tool and its implementation.

2.1. Gpower-normal Distribution

With the goal mentioned above, a new distribution family named gpower-normal is presented and its main properties are stated. We are considering data that, when transformed by gpower, become normally distributed.

Definition 1.

A random variable Y has a gpower-normal distribution if for some $\mu \in \mathbb{R}$, $\sigma >0$, $p\in \mathbb{R}$ its density function is given by

$$f_Y\left(y;\mu ,{\sigma }^2,p\right)=\frac{1}{K\sqrt{2\pi {\sigma }^2}}\frac{{\left(y+\sqrt{y^2+1}\right)}^p}{\sqrt{y^2+1}}exp\left(-\frac{1}{2{\sigma }^2}{d}_p\left(y,\mu \right)\right),y\in \mathbb{R}$$

(2)

where the normalizing constant $K=1-\Phi \left(\frac{-1/p-gpower\left(\mu ,p\right)}{\sigma }\right)ifp>0$, $K=1ifp=0$ and $K=\Phi \left(\frac{-1/p-gpower\left(\mu ,p\right)}{\sigma }\right)ifp<0$. Φ is the cumulative distribution function of a standard normal variable and $d_p\left(y,\mu \right)={\left(gpower(y;p)-gpower(\mu ,p)\right)}^2$ is the deviance (see Section 2.2). Density 2 will be denoted as $GPN(\mu ,\sigma ,p)$, where GPN stands for gpower-normal distribution.

The next theorem gives the main property of gpower-normal variables: after a gpower transformation they become truncated normals (TN). Recall that if X is a TN variable, its density is given by $f_X\left(x,{\mu }_X,{\sigma }^2\right)=\frac{1}{K}\frac{1}{\sqrt{2\pi }\sigma }exp\left(-\frac{1}{2}{\left(\frac{x-{\mu }_X}{\sigma }\right)}^2\right)I_{\left(a,b\right)}$, where I is the indicator function, $\left(a,b\right)=\left(-1/p,\infty \right)\mathrm{if}p>0$, $\left(a,b\right)=\left(-\infty ,\infty \right)\mathrm{if}p=0$ and $\left(a,b\right)=\left(-\infty ,-1/p\right)\mathrm{if}p<0$; ${\mu }_X\in \mathbb{R}$ and we will denote that $X\sim TN\left({\mu }_X,{\sigma }^2,a,b\right)$ (see Dhrymes [13]).

Theorem 1.

Let $Y\sim GPN(\mu ,\sigma ,p)$, then the transformed variable $X=gpower\left(Y;p\right)$ has a truncated normal distribution $TN\left({\mu }_X,{\sigma }^2,-1/p,\infty \right)$ if $p>0$, $TN\left({\mu }_X,{\sigma }^2,-\infty ,-1/p\right)$ if $p<0$ and normal distribution $N\left({\mu }_X,{\sigma }^2\right)$ if $p=0$ with ${\mu }_X=gpower(\mu ,p)$.

The proof of this theorem can be seen in Appendix A. Figure 1 shows the flexibility of this distribution family across different values of the parameters, including symmetric and heavy tailed densities.

2.2. Relationship between Gpower-normal Models and Pseudo-dispersion Models

Gpower-normal models are a special case of a general family of distributions called pseudo-dispersion models defined by [14]. It is proven in Appendix B that the densities defined by (2) belong to this family. Expressions for their deviance and unit variance functions are also obtained. From a theoretical point of view, this is interesting because there are very few examples of pseudo-dispersion models known in the specialized literature.

2.3. Quantiles and Moments

A straightforward method to obtain quantiles such as the median and quartiles is considered. These quantiles will be useful in the graphical examination of the model fit to a data set, through quantile–quantile plots.

Let Y be a gpower-normal random variable and $X=gpower\left(Y;p\right)$ the transformed variable distributed as $TN\left({\mu }_X,{\sigma }^2,-1/p,\infty \right)$. Let $x_{\alpha }$ be the α-quantile for X, to obtain an expression for the quantiles of a truncated normal distribution we proceeded as follows:

$$\begin{array}{cc}\hfill P\left(X\le x_{\alpha }\right)& =P\left(X^0\le x_{\alpha }\right)-F_{X^0}\left(-1/p\right)\hfill \\ & =\frac{1}{K}\left(\Phi \left(\frac{x_{\alpha }-{\mu }_X}{\sigma }\right)-\Phi \left(\frac{-1/p-{\mu }_X}{\sigma }\right)\right)\hfill \\ & =\frac{1}{K}\left(\Phi \left(\frac{x_{\alpha }-{\mu }_X}{\sigma }\right)-\left(1-K\right)\right)=\alpha ,\hfill \end{array}$$

where $X^0$∼$N\left({\mu }_X,\sigma \right)$ is the corresponding normal variable with cumulative distribution $F_{X^0}$ and Φ is the cumulative distribution of the standard normal. Now, clearing up

$$x_{\alpha }=\sigma {\Phi }^{-1}\left(K(\alpha -1)+1\right)+{\mu }_X.$$

Also

$$\alpha =P\left(X\le x_{\alpha }\right)=P\left(gpower\left(Y\right)\le x_{\alpha }\right)=P\left(Y\le gpowe{r}^{-1}\left(x_{\alpha }\right)\right)$$

then the α-quantile for Y is

$$y_{\alpha }=gpowe{r}^{-1}\left(x_{\alpha }\right)$$

(3)

and its value can be obtained from the standard normal distribution. This procedure will be applied in Section 3.2.

Expressions for the moments of a gpower-normal family can be expressed in terms of the truncated normal density function (see Appendix C). However, these expressions are not easy to handle, and their convergence for different values of the parameters remains to be analysed. In order to avoid the difficulties in moment calculations (e.g., mean and variance) for models given in Definition 1, we propose the alternative use of quantiles.

2.4. Parameter Estimation

The gpower-normal models have three parameters to be estimated. They are related to the corresponding $TN$ model parameters as it has been detailed in Section 2.1. We propose a combined profile likelihood and maximum likelihood approach to estimate the parameters. The five steps of the proposed estimation approach are described below:

• Given a data set represented by vector $\mathbf{y}=(y_1$, $y_2$, …, $y_n)$, to obtain a profile likelihood for the power p, we consider a grid of values $p_0$, $p_1$, …, $p_k$.
• For each $p_j$, $1\le j\le k$ the transformed data ${\mathbf{x}}_{p_j}$ are calculated as ${\mathbf{x}}_{p_j}=gpower\left(\mathbf{y},p_j\right).$
• Then, for each $p_j$, the corresponding ${\mu }_{p_j}$ and ${\sigma }_{p_j}$ are estimated, maximizing the likelihood function of the truncated normal variable.
• Then, $p_j$, ${\mu }_{p_j}$ and ${\sigma }_{p_j}$ are used to obtain the log-likelihood function of $\mathbf{y}$ whose density was given by (2):

  $$ln{f}_{\mathbf{y}}\left(\mathbf{y};\mu ,{\sigma }^2,p\right)=nln\frac{1}{K\sqrt{2\pi {\sigma }_{p_j}^2}}+{\sum }_{i=1}^n\left(lnp\left(y_i+\sqrt{y_i^2+1}\right)-ln\sqrt{y_i^2+1}\right)-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^n{d}_p\left(y_i,{\mu }_{p_j}\right).$$
• Finally, p is chosen as the one that maximizes the log-likelihood in the grid: $\widehat{p}=\underset{1\le j\le k}{max}ln{f}_{\mathbf{y}}\left(\mathbf{y};{\mu }_{p_j},{\sigma }_{p_j}^2,p_j\right)$

The method described above is applied in Section 3.2 and an implementation has been written in R language [15]. The codes are available from the authors upon request.

3. Results

To highlight the potential of the distribution family to model data with different skewness and kurtosis, we present several density functions and their quantiles. Also, the model fit is illustrated with real data sets.

3.1. Some Examples of Gpower-normal Densities

By way of example, Figure 1 shows the flexibility of this distribution family across different values of the parameters, including symmetric and heavy tailed densities. Graphic representations of these densities are exhibited for two different values of p (0.05 and 0.5) and three different values of σ (1, 5 and 10); always with $\mu =0$.

It can be seen that kurtosis decreases as σ grows up and asymmetry grows with p.

Table 1. Some α quantile values for different choices of the parameters in Expression 2: $p=0.05,0.50$, $\sigma =1,5,10$ always with $\mu =0$.

		p = 0.05			p = 0.50
α	σ = 1	σ = 5	σ = 10	σ = 1	σ = 5	σ = 10
0.01	−11.78	−3.65 × 7	−1.63 × 22	−166.37	−508.68	−182.81
0.05	−5.39	−3.98 × 104	−2.73 × 1011	−13.84	−20.83	−7.25
0.10	−3.49	−2.26 × 103	−4.64 × 107	−5.67	−5.17	−1.33
0.50	0.00	3.97 × 10−4	0.57	0.06	4.25	14.16
0.90	3.17	259.00	2.17 × 104	2.34	22.75	78.52
0.95	4.65	981.00	1.74 × 105	3.04	31.37	109.85
0.99	8.92	9590.00	5.24 × 106	4.48	51.85	185.57

contains some quantile values corresponding to the distributions displayed in Figure 1 confirming the observations made above for those displays.

3.2. Real Data Applications

As it was mentioned in Section 1, the proposed density family was originally motivated by the modelling of microarray intensities, but its application is more general. Here, we present some examples: the first one corresponds to intensities of microarray data, the second one to concentrations of ammonia and the third one to magnetic contamination indices. For these examples, the parameter estimators are obtained by the method described in Section 2.4. Then, the data fit to the corresponding estimated model is shown graphically in three ways: (1) by the overlap of the data histograms with the density curves, (2) by the overlap of the empirical distribution curve with the adjusted model cumulative distribution function, and (3) by quantile–quantile (Q–Q) plots, that display the ordered data in the horizontal axis and the corresponding quantiles of the estimated model distribution.

It is also evaluated applying Kolmogorov–Smirnov tests to compare the estimated gpower-normal and the glog-normal distributions. The resulting p-values should be taken as descriptive comparisons of the two distributions and as a complement of the Q–Q plots.

Example 1.

The first set of data corresponds to 30 intensities of one gene and they can be seen in Appendix D. These data were selected from the Yale University MAQC project and downloaded from [16].

The values for p, μ, and σ were estimated by the profile likelihood method (Section 2.4), resulting $\widehat{\mathrm{p}}$ = 0.23, $\widehat{\mu }$ = 1052 and $\widehat{\sigma }$ = 3.46. From these estimations, the 0.25, 0.50 and 0.75 quantiles were obtained according to expression (3) resulting in 690, 1052 and 1545, respectively. As can be seen in Figure 2, this set of data is well fitted by the gpower-normal model with the parameters given above. It is confirmed by the Kolmogorov–Smirnov test, with a p-value = 0.7057 for the gpower-normal model and 0.4934 for the glog-normal model.

Example 2.

Two sets of data representing concentrations of ammonia in the UK for the years 2005 and 2011 were downloaded from [17] and they are given in Appendix D. Applying again the methodology detailed in Section 2.4, we obtained, for the year 2005, an optimum value of 0.20 for the parameter p. The corresponding, μ, σ and the 0.25, 0.50 and 0.75 quantiles were 2.11, 0.70 , 1.39, 2.11 and 3.04 respectively. For the year 2011, the optimum p was 0.16 and the corresponding $\widehat{\mu }$, $\widehat{\sigma }$ and the 0.25, 0.50 and 0.75 quantiles were 2.93, 0.68, 2.00, 2.93 and 4.15 respectively. Observe that the 2011 interquartile interval shows a level and amplitude increase with respect to 2005.

As in the previous example, it can be seen in Figure 3 and Figure 4 that the fitted gpower-normal models with the parameters given above fit them quite well.

In the Q–Q plot for the year 2011 data, the highest three observations are higher than the expected model quantiles, suggesting an increase in ammonia concentration. The Kolmogorov–Smirnov test confirms the goodness-of-fit with a significance of 0.3411 for 2005 and 0.4993 for 2011. For the corresponding glog-normal model, p-values of 0.1168 and 0.2327 were obtained.

Example 3.

When measuring magnetic contamination, the index can take an asymmetric positive distribution, with some negative values observed because of air presence. In this example, data from Mar del Plata, Argentina (see [18]), are fitted with a gpower-normal density. Data can be seen in Appendix D. Applying again the methodology detailed in Section 2.4, we obtained an optimum value of 0.19 for the parameter p. The corresponding $\widehat{\mu }$, $\widehat{\sigma }$ and the 0.25, 0.50 and 0.75 quantiles were 8.05, 1.64 , 3.98, 8.05 and 14.91 respectively.

As in the previous examples, it can be seen in Figure 5 that the fitted gpower-normal model with the parameters given above fits them quite well. The fit was evaluated applying a Kolmogorov–Smirnov test (p-value = 0.57 for the gpower-normal model and 0.29 for the glog-normal model).

4. Discussion

A new family of distributions named gpower-normal, indexed by $p\in \mathbb{R}$, has been defined. Variables whose distribution belongs to this family become normal or truncated normal when a gpower transformation is applied. From a practical point of view, the truncation is very often negligible, such is the case of data in Example 1. The estimated parameters were $\widehat{p}=0.23$ , $\widehat{\mu }=1052$ and $\widehat{\sigma }=3.516$ and thus the truncation constant given in Definition 1 is almost 1, meaning that no truncation is necessary in this case. Similar results were obtained for Examples 2 and 3.

The gpower-normal family can model data that include non-positive values, as required for microarray analysis.

To estimate the model parameters, a combined maximum likelihood method is proposed and it is successfully applied to real data. It enables direct calculations of quantiles for which simple expressions are given. Thus, position and scale measures (i.e., the median and the interquartile range) can be easily obtained, overcoming the difficulties found in moment calculations.

This allows the use of the estimated distribution medians as summary measures for gene intensities. Therefore, to compare the gene intensities in two biological conditions, it seems adequate to compare the corresponding medians of the estimated distributions.

This paper extends the results previously presented by Leiva and coauthors [11] and therefore the new distribution family offers a larger set of models. Considering the Leiva et al. proposal as a standard alternative, the new family can fit data for which their proposal might be not flexible enough. As expected, in all examples, the Kolmogorov–Smirnov tests indicated a closer fit of the data to the gpower-normal distribution in comparison with the glog-normal fit.

In addition, it has been proven that the gpower-normal family is a special case of pseudo-dispersion models, inheriting all the good properties of these models, such as asymptotic normality for small variances. It should be pointed out that very few examples of pseudo-dispersion models have been reported in the literature. The obtained p-values should be considered as descriptive values only and as a complement to the Q–Q plot, which provides a visual impression of the fit.

5. Conclusions

We have presented a family of distributions to model asymmetric data with positive and negative values. It allows investigators to model microarray data, preserving their original scale. Although the construction of this family of distributions was motivated by microarray data, it is important to notice that the model is also useful in other environments, such as chemical and magnetical contamination data. A method has been given to estimate quantiles and in our examples, specific quantile–quantile plots were effective tools for a visual evaluation of the data fit to the estimated models.