{'fraction_non_alphanumeric': 0.02980120372059092, 'fraction_numerical': 0.007514134597847893, 'mean_word_length': 4.8431372549019605, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Bayesian model for microarray datasets merging

The aggregation of microarray datasets originating from different studies is still a difficult open problem. Currently, best results are generally obtained by the so-called meta-analysis approach, which aggregates results from individual datasets, instead of analyzing aggre-gated datasets. In order to tackle such aggregation problems, it is necessary to correct for interstudy variability prior to aggregation. The goal of this paper is to present a new approach for microarray datasets merging, based upon explicit modeling of interstudy variability and gene variability. We develop and demonstrate a new algorithm for microarray datasets merging. The underlying model assumes normally distributed intrinsic gene expressions, distorted by a study-dependent nonlinear transformation, and study dependent (normally distributed) observation noise. The algorithm addresses both parameter estimation (the parameters being gene expression means and variances, observation noise variances and the nonlinear transformations) and data adjustment, and yields as a result adjusted datasets suitable for aggregation. The method is validated on two case studies. The first one concerns E. Coli expression data, artificially distorted by given nonlinear transformations and additive observation noise. The proposed method is able to correct for the distortion, and yields adjusted datasets from which the relevant biological effects can be recovered, as shown by a standard differential analysis. The second case study concerns the aggregation of two real prostate cancer datasets. After adjustment using the proposed algorithm, a differential analysis performed on adjusted datasets yields a larger number of differentially expressed genes (between control and tumor data). The proposed method has been implemented using the statistical software R 1, and Bioconductor packages 2. The source code (valid for merging two datasets), as well as the datasets used for the validation, and some complementary results, are made available on the web site


Introduction
Microarray technology provides high throughput measurements of messenger RNA levels of thousands of genes in tissue samples. Since its introduction almost twenty years ago, it has found applications in many aspects of molecular genetics and functional genomics, from the discovery of basic biological mechanisms to the classification of diseases into subgroups and the prediction of disease outcome.
However, although important progress has been made, that has turned lab experiments into industrially standardized protocols, the technology is still facing significant reproducibility problems (see e.g. the introduction of [5] for a more detailed account). Therefore, comparing results from different studies is generally difficult, and extreme care is needed if one wishes to merge different datasets into a single one to improve reliability and generalizability of the results.
Several microarray study aggregation techniques have been proposed. Among them, two approaches have received significant attention. The first one (see for example [4]) proposes to replace numerical expression data with ranked data, in order to correct for study-dependent normalization effects. This implicitely assumes that the ranks are essentially respected in the considered datasets; unfortunately, such an assumption fails to be true in many situations, for example when a large number of weakly expressed genes are present, the fluctuations of which induce large deviations in ranks. The second classical approach, termed Meta-Analysis, aggregates the results of statistical analyzes from the different datasets rather than aggregating the datasets themselves (see for example [1] for a review). It has been argued by various authors (see for example [1,2]) that such meta analysis approaches generally outperform classical statistical analysis directly performed on aggregated datasets.
We propose in this paper a new approach for pre-processing gene expression datasets prior to aggregation, as an alternative to meta-analysis based techniques. The idea is to start with an explicit model for gene expression data, and an explicit model for the distortions induced by the different experiments to be merged. In a Bayesian framework, those two models are used as priors. Assuming normally distributed additive observation noises leads to a closed form expression for the posterior probability.
This model is simple enough to yield a fairly standard estimation algorithm, exploiting elementary linearization and optimization techniques. The latter yields estimates for the model parameters, together with adjusted gene expression values in which both the non-linear distortions and the variances inhomogenities have been corrected, therefore suitable for aggregation and further statistical analysis.
The model and the algorithm are validated in both synthetic datasets (i.e. real datasets with an artificial distortion and additional noise), for which the method is shown to be able to correct for the distortion and noise (to some extent). We then apply the approach to real prostate cancer expression datasets, for which two original datasets are aggregated after correction.
Preliminary results of this work have been given in conference proceedings [8].

Approach
The approach we develop here is based upon an explicit modeling of the gene expressions, and of the distortion introduced by the several studies. For the sake of simplicity we limit ourselves to a simple model for logarithms of gene expressions, assuming independent normally distributed values (with unknown means and variances). We also assume that each study gives rise to a non-linear distortion (called an observation function) of gene expressions, and an additive zero-mean Gaussian white noise, of unknown variance. Our approach mainly relies on two steps. First, estimate the parameters of the model (including the nonlinear observation functions), then estimate adjusted gene expression values, that can be aggregated into a single dataset.
More precisely, we rely on the following model for logarithms of measured gene expression data. Assume we are given several datasets k = 1, . . . K, involving g = 1, . . . N g genes (if necessary, after restriction to a suitable subset of common available genes). For each dataset k, several arrays (hereafter termed conditions c = 1, . . . N k c are available. We denote by y k gc the logarithms of measured gene expression levels, which we shall call throughout this paper the observed gene expressions. The ingredients of the model are the following: • The intrinsic gene expression values x = {x k gc } are assumed to be independent realizations of i.i.d. normally distributed random variables X g ∼ N (µ g , σ 2 g ) (whose distributions do not depend on the study k). We write x k gc = µ g + δ k gc .
• The observed gene expression values y = {y k gc } are obtained by a study-dependent nonlinear transformation f k of the intrinsic expression values, with additive Gaussian noise, Y k g = f k (X g )+ U k , with U k ∼ N (0, τ 2 k ). Componentwise, the observations therefore read • The observation functions : f = {f k , k = 1, . . . K} are assumed to be monotonic, smooth functions ( |f ′′ k (x)| 2 dx < ∞), and are given a spline-type prior distribution Conditionally to the observation functions f k and the intrinsic gene expressions x k gc , the observed gene expressions y k gc are therefore independent and normally distributed, with means f k (x k gc ) and variances τ 2 k . Applying the Bayes rule, it is easily seen that the log posterior probability distribution is (up to an additive constant) equal to the sum of three components The approach we develop here aims at maximizing numerically the above log-posterior probability, estimating the gene average expressions µ g and variances σ 2 g , the observation noise variances τ 2 k and the observation functions f k , and solve the regression problem for estimating intrinsic gene expression values x k gc .

Methods
We describe in this section the estimation method used in our approach, together with the corresponding algorithm. The first step of the algorithm is the estimation of the model parameters: means, variances and observation function. The chosen approach is an iterative one, each parameter is estimated at the time. After the estimation, adjusted gene expression datasets are available, that can be aggregated and analyzed. and decouples further into the following K independent optimization problems: for k = 1, . . . K, solve

Estimation
The latter are actually smoothing spline regression problems (see e.g. [11] for a detailed account), for which efficient algorithms are available and implemented in standard software packages. It is worth noticing that spline functions being piecewise polynomials, once the spline has been estimated, its derivative is readily available.
Remark 1 These estimations are performed on a set of genes with small variance across samples in each experiment, termed below invariant gene set. This gene set must also be representative of the range of values of the average gene expressions. Practically, the invariant gene set is determined as follows: • Split the range of values of the observed gene expressions into intervals.
• Within each interval, select a given percentage of genes with smallest variance.

Mean and variances estimation
The average gene expressions µ g are re-estimated at each step of the algorithm as weighted sample averages of the estimated gene expressions: at iteration t, where N c = k N k c is the total number of arrays. The variance component estimation is a difficult task here, as many gene variances σ 2 g are to be estimated. The I-MINQUE (or REML) approach (see e.g. [6]) is a natural choice, which is unfortunately not suitable here as it produces negative estimates for the gene variances. To estimate the latter, we choose sample estimates from the initialization (i.e. pre-adjusted expression data, see Section 3.1.4 below). To estimate the observation noise variances τ 2 k , we resort to the I-MINQUE approach, applied to the invariant gene set (defined in Remark 1).

Remark 2
Our simulations (see the E.Coli example below) show that when the observation noise is large, the corresponding variances τ 2 k can be underestimated by our procedure, which in turn results in an overestimation of the gene variances. To overcome this problem, a possible solution is to regularize the gene estimates by adding a positive constant λ 2 to the estimated noise variances: modified estimates of the formτ 2 k = τ 2 k + λ 2 are considered. In practical situations, since the observation noise is unknown, the regularization parameter λ can be used as a tuning parameter.

Adjusted gene expression values estimation
Suppose now that the observation functions f k are known, the minimization of the log posterior probability in (1) leads to the optimization of which splits into a sum of N g decoupled terms: the estimation of the intrinsic gene expression values x k gc reduces to minimizing for each g = 1, . . . N g the quantity Due to the non-linearity of the observation functions f k , no closed-form expression exist for the solution, and we resort to an iterative numerical algorithm. We assume that the mean µ g and variances σ 2 g and τ 2 k are known, as well as the observation functions f k . At iteration t, assume that the previous estimate x k gc (t−1) of the gene expression values is available. A linearization of the observation functions f k in the neighborhood of x k gc (t − 1) yields the first order approximation x k gc (t) = x k gc (t − 1) + ǫ k gc , and The update x k gc (t) can therefore be obtained by optimizing with respect to ǫ k gc the above expression, which yields Then 0 ≤ α k gc ≤ 1, the bounds being attained in the extreme cases (no noise, or constant f k ). This yields the update rule x k i.e. a weighted average of the mean µ g and the contribution of observations.
Remark 3 This is similar to empirical Bayes type update rules, the difference being that the weights depend upon the observations, due to the nonlinearity of observation functions.
The balance between the contributions of the mean and the observations is controlled by the signal to noise ratio (SNR) σ 2 g /τ 2 k . The larger the SNR, the smaller the α parameter and the larger the contribution of the observations: adjusted values are closer to the observed ones. Conversely, for small SNR values, adjusted gene expressions are closer to the means µ g .

Remark 4
It may be shown that this update rule is actually equivalent to a variable step gradient descent method: , where the stepsize depends on the gene g, the condition c, the study k and the iteration index

Initialization.
The estimation procedures described above are used within an iterative algorithm. The latter needs either initial values for the intrinsic gene expressions, or the observation functions, since only data y is available. We use an estimate for the reciprocal of the observation functions (the so-called rectification functions) provided by the approach described in [7]. The estimation of rectification function is formulated in [7] as another smoothing spline problem: optimize, with respect to the mean gene expressions µ g and the rectification functions ϕ k the quantity The problem is solved by an iterative algorithm, in the same spirit as the approach described here. Besides the fact that the approach of [7] estimates rectification functions instead of observation functions, the main difference with the current approach is that the latter models observation noise, which is not the case of the former. The approach developed in the present work is more complex, but is able to correct for study-dependent observation noise.

Algorithm
The proposed approach can be summarized as follows: • Initialization: Start from a first estimate x k gc (0) for the intrinsic expression values, provided by rectification function estimates (see section 3.1.4). Estimate the gene average expressions µ g and variances σ 2 g as in section 3.1.2, and the observation functions as in section 3.1.1.
• Iteration t: estimates x k gc (t − 1) are available. * Re-estimate the gene expressions x k gc (t) as in 3.1.3. * Update the mean gene expressions µ g (t) as in 3.1.2.
The output of the algorithm consists in estimatesx = {x k gc } for the expression datasets x = {x k gc }, to be exploited for further analyses, together with estimates for the means µ g , variances σ 2 g and τ 2 k , and the observation functions f k .
The algorithm was implemented using the R statistical environment, from which we used the smoothing spline function smooth.spline. Bioinformatics related functions from the Bioconductor package were also used.

Validation: differential analysis
Even though differential analysis is not the only application we have in mind, it provides a convenient setup for validating the proposed approach. In the numerical results we discuss below, differentially expressed genes were searched for in original (distorted) datasets and adjusted datasets, using t-test with FDR correction for multiple testing. We used functionalities from the Bioconductor package, namely the MTP function from the multtest package.

Results and discussion
The proposed approach has been tested and validated on several problems, of increasing complexity. We first verified that the algorithm performs well on articifial data, i.e. random data simulated according to the above model. The results are quite satisfactory, but since there's not much to discuss about such simulations we refrain from reporting on them here, and prefer to focus on results closer to practical situations. We report below on two different validations of the approach. The first one concerns semi-artificial data, namely real data, with clear biological outcome, which are artificially distorted according to the distortion model above. We show that in such a situation, the proposed approach allows one to almost completely correct for the effect of the distortion, and recover the biological meaning of the data. The second validation is performed on two real human gene expression datasets, focusing on prostate cancer.

Real data with artificial distortions
A test was performed using artificial observation functions, applied to real data. Namely, we chose a dataset with well understood biological outcome, splitted it into two well balanced subsets denoted by E 1o and E 2o ("o" stands for "original", see below for details) and applied to the two subsets two different observation functions, before adding Gaussian observation noise (yielding distorted subsets E 1d and E 2d . After adjustment, datasets denoted by E 1a and E 2a were obtained, that were aggregated. The goal was to study the impact of the deformation induced by the observation functions and the noise (which were chosen so as to hide the biological effects), and the ability of the algorithm to perform a sensible correction. E. Coli expression data from [3] were used. The data include expressions of 7295 genes under two different situations (20 aerobic and 22 anaerobic). They exhibit a clear variability between the two biological situations, and are thus particularly interesting since they provide a good test of the ability of our approach to recover such a variability after distortion and correction.
Two subsets were created, both involving randomly chosen 10 aerobic and 11 anaerobic conditions. Non-linear transformations f 1 and f 2 were applied to the two so-created subsets (after standardization), and observation noise with variances τ 2 1 and τ 2 2 was added. We report here on the particular case    Table 1 for the corresponding observation noise variances), the sign indicates over or under expression in aerobic conditions. Column 2: numbers of differentially expressed genes in the original dataset. Columns 3 to 7: numbers of differentially expressed genes in distorted (D) and adjusted (A) datasets, and in the intersections E 1o ∩ E 2o , E 1d ∩ E 2d and E 1a ∩ E 2a that are present in the original dataset (the numbers between parentheses indicate the total numbers of differential genes).
with various values of the tuning parameter α. Actual variance values, together with the corresponding estimated values, are displayed in Table 1, together with the estimated average gene variance, denoted by σ 2 g . As can be seen, the estimates are fairly good, except for the last situation, which corresponds to the value α = 20, presumably too large for our procedure (results continue to degrade for higher values of α, which probably originates from the degradation of the gene variance estimates).
To illustrate graphically the effect of the adjustment, we display in Fig. 2 the first factorial plane obtained by (normalized) principal component analysis with original data, distorted data, pre-adjusted data (using the method of [7], used here as initialization, see Section 3.1.4), and adjusted data. The projection of the output datax k gc onto the first factorial plane turns out to reproduce fairly accurately the projection of the original data x k gc (before distortion). The processing has therefore permitted to recover the biological features as the first source of variability. This is confirmed by visual inspection of the projections on other factorial planes (not shown here).
A more quantitative validation can be obtained through differential analysis. We searched for genes differentially expressed between aerobic and anaerobic conditions, on initial data, distorted data and adjusted data, for various values of observation noise. For the sake of comparison with meta-analysis based approaches, we also performed differential analysis on the two subsets E 1 and E 2 (in original, distorted and adjusted situations) and computed the intersection of the two estimated gene subsets in each situation. The goal of the comparison is to asses the ability of each approach to recover original differentially expressed genes.
The results are displayed in Table 2, which gives the numbers of differentially expressed genes on the original dataset (first column), and the number of original differential genes that are recovered on the distorted dataset (second column) and the adjusted dataset (third column); total differential gene numbers are given between parenthesis.
A first conclusion is that the distortion significantly degrades the differential analysis, the degradation being obviously bigger when the observation noise is higher. The processing improves the situation, which is reflected by an average increase of the number of the recovered original differential genes ranging between 11% and 20% (depending on the noise level).
For the sake of comparison, we also display in the last three columns numbers of differential genes obtained by intersecting the results of individual differential analysis on the two sub-datasets, in the three considered situations. It is worth noticing that the loss of differential genes is in this case really huge, which is for us a strong motivation for aggregating datasets rather than intersecting differential analysis results.

Real data
The proposed approach was also tested on the two publicly available datasets of prostate cancer expression data reported in [9] and [10]. Both datasets originate from Affymetrix HG-U95av2 arrays.
We started from raw expression data, and normalized them separately using the standard RMA normalization procedure from the Bioconductor package. After normalization, the dataset of [9] appeared to exhibit much larger variability than the [10] dataset. We thus performed the following pre-processing: we extracted a subset of arrays whose correlations to the median array exceed 90%, and reduction to common genes). After pre-processing, the two-datasets consist in respectively 61 (32 tumor and 29 normal) and 86 (37 tumor and 49 normal) conditions, with 12625 genes, from which we extracted balanced sub-datasets (29 tumor and 29 normal for each)  The proposed algorithm was run on these two datasets. The result of the processing is shown in Fig. 3 (estimates for the observation functions) and in Fig. 4 (estimated variances, for the observed and adjusted data). We first notice that the estimated observation functions are in this case fairly regular, no large distortion has been found. It can be seen on Fig. 3 that the Singh data (termed E 1 ) have been shrunk more importantly than the Stuart data (E 2 ). This can be interpreted in terms of the observation noise. Indeed, the observation noise was found to be significantly larger in the Singh Finally, it is clear from Fig. 4 that the proposed approach has done quite a good job for adjusting the gene variances: variances in adjusted Singh and Stuart datasets are coherent. The interstudy variability has therefore been considerably reduced.
As in the previous case study, to assess more quantitatively the performances of the proposed approach, differential analysis was performed on the real dataset, and the adjusted dataset. After filtering out the 30% least variable genes, we used t-test based decision rule, with FDR correction for multiple testing (α = 5%, 5000 bootstrap samples). Differentially expressed genes were seeked for the real dataset and the processed dataset, as well as the individual subsets E 1 (Singh) and E 2 (Stuart), and the corresponding adjusted subsets E 1a and E 2a . The results are summarized in Table 3  As before, a striking first result is the fact that the intersection of the differential genes found from E 1 and E 2 individually is extremely small, as well as the the intersection of the differential genes found from E 1a and E 2a individually. Quite a large number of differential genes appear to be lost by this procedure. Besides, as was also seen in the above E. Coli simulation study, the adjusted datasets yield a larger number of differential genes than the original datasets. A closer examination of the results of the differential analysis leads to the following conclusions:  Figure 5: P-values of differential genes for the original and the adjusted datasets: Positively differential genes in the original (O+) and adjusted (A+) datasets, and negativeley differential genes in the original (O-) and adjusted (A-) datasets differential genes. Even if one focuses on differential genes with small p-values (say, less than 1%), the adjusted dataset yields a significant number of "new" genes, including for example MAP-KAPK3, POLR2H, BDH1, PYCR1, PDIA4, SLC25A6, ZMPSTE21, ENTPD6 for the positively differential genes (i.e. overexpressed in tumors), and ADD1, VAT1, PPP2CB, RBPMS, FSCN1, ARHGEF4, GRK5, FEZ1, HTRA1, RBMS1, CETN2 and OSR2 for the negatively differential genes.
A precise analysis of the role of these genes in prostate cancer goes far beyond the scope of the present work. Let us simply notice that most of these genes have already been reported in the prostate cancer literature. This makes them potentially interesting.
• Globally, the p-values of differential genes are smaller in the adjusted dataset than in the original dataset. This appears clearly in Fig. 5 where they are displayed in boxplot forms.
• Conversely, less genes were found differential in the original dataset only; in addition, the corresponding p-values are generally close to the critical value chosen in this study (5%). A closer examination shows that their variance in the adjusted dataset is smaller than in the original one. According to Remark 3, this presumably originates from a low signal to noise ratio σ 2 g /τ 2 k , or a large value of the observation function derivative f ′ (µ g ) (in the present approach, these quantities govern the parameter α k gc which controls the variance in the adjusted dataset).

Conclusion
We have proposed in this paper a new approach for aggregating gene expression datasets originating from different studies. The rationale of the method is to adjust the datasets prior to merging, using a Bayesian modeling. The model we develop is specifically adapted to microarray data, and assumes that the inter-study variability originates from non linear distortions, combined with study dependent observation noise. The inter-study variability model is combined with a Gaussian model for the gene expression values (after logarithmic transformation) in a Bayesian framework. Our results on artificial and real data show that the method is sound, and capable of correcting for study-dependent distortion and observation noise. For the sake of simplicity, we focused on the case of two datasets to be merged, however the model is versatile enough to accomodate arbitrary numbers of datasets. The results on real data discussed here are quite encouraging, however the approach will have to be tested further on a larger scale to be validated more thoroughly.
Several aspects are to be investigated further. Among them, the problem of variance components estimation is an important one, as the estimated gene variances have a key impact on the quality of the results. For the sake of simplicity we have limited ourselves here to a simple procedure, involving sample estimates of gene variances from pre-adjusted data. More sophisticated approaches will have to be investigated.