A Flexible Microarray Data Simulation Model

Microarray technology allows monitoring of gene expression profiling at the genome level. This is useful in order to search for genes involved in a disease. The performances of the methods used to select interesting genes are most often judged after other analyzes (qPCR validation, search in databases...), which are also subject to error. A good evaluation of gene selection methods is possible with data whose characteristics are known, that is to say, synthetic data. We propose a model to simulate microarray data with similar characteristics to the data commonly produced by current platforms. The parameters used in this model are described to allow the user to generate data with varying characteristics. In order to show the flexibility of the proposed model, a commented example is given and illustrated. An R package is available for immediate use.


Introduction
Microarray is now a mature technology in the field of molecular biology for monitoring of gene expression profiling at the genome level. Using two-color microarray technology, the biological activities of two samples are compared on an array on which thousands of specific deoxyribonucleic (DNA) sequences are printed or synthesized in situ. Two fluorescent dyes (generally red and green) are used for labeling of the samples for hybridization. Then, after washing, scanning and quantification of images corresponding to the red and green fluorescence of the array, numerical values, or intensities are associated with each probe (gene). These values are normalized to correct for undesirable technical variations [1,2]. Using one-color microarray technology, each sample is hybridized on one array. One fluorescence color is used, and expression intensities are obtained for genes on the arrays used. These values should also be normalized [3]. For a given gene, a ratio of values from two samples is used as a measurement of its expression change. Throughout this paper, we use logarithm two scale to transform the values. Hence, we use log 2 intensities and log 2 ratios for sample intensities and their ratios. In general, log 2 ratios are used for two-color microarray data, while log 2 intensities are used for one-color microarray data.
In a screening study, the goal is to select genes with differential expression from data whose characteristics are unknown. To study the performance of gene selection methods, we often use synthetic data, which can also serve in developing new methods. To properly play their role, synthetic data must resemble as closely as possible the real data they represent. This is achieved by simulating the physical phenomena through a model. The parameters of this model are approximations of the physical laws that govern the observed phenomenon, and a good knowledge of the physical laws is therefore required to obtain a good model. A complex model that takes into account several components of a phenomenon may become less flexible and may be difficult to modify in order to include new factors. In contrast, a simple model can be effective and easy to modify to take into account an unexpected situation.
One way to generate synthetic data is to use real microarray values as seeds [4][5][6][7]. In [4], data were generated using a normal distribution and a real microarray dataset. This dataset was used to estimate some hyperparameters, and the percentage of the differentially expressed (DE) genes was fixed; see details in [4]. The simulation data obtained in [4] may be useful for statistical methods, but they differ from observed data, since the log 2 -intensities obtained vary between the unrealistic bounds 16 and 30. In [5], two samples are selected from the control and the test samples. The measurements associated with their genes are modified (exchange of values between the two samples. . . ) in order to obtain two statistically undistinguishable samples. Finally, a given number of down-and up-regulated genes is used in this dataset. This procedure necessitates a real microarray dataset, which cannot be available. The procedure proposed in [6] is a modular system, including all steps of the microarray technology. This includes slide layout, hybridization, scanning and image processing. Many models already available in the literature are used with new ones in [6]. The system of [6] results in data as close as possible to real biological data, but is quite complex, and the large number of its parameter settings may discourage its use. In [7], two simulation methods are used. In the first, the number of DE genes and the number of levels of changes for these genes are fixed. Then, for each gene, a mean and a standard deviation are drawn from a uniform distribution. Finally, these data (mean and standard deviation) are used as normal distribution parameters to get expression values for the genes. In the second simulation method, the mean and standard deviation of test (µ t , σ t ) and control (µ c , σ c ) samples are estimated from observed real data and are used as normal distribution parameters to get expression values for all genes. There is no flexibility in the number of the DE genes of the method in [7]. A hierarchical model is used in [8,9], where the variance of each gene is simulated using two parameter settings and the χ 2 distribution. This variance is then used to generate values from a normal distribution. DE genes are finally defined using levels compared to a threshold. The procedure used in [8,9] allows one to generate data with parameters derived from assumptions about real data. However, the percentage of the DE genes depends on a parameter that is difficult to control a priori. For some of the above methods, there is no distinction between the number of weakly expressed genes and those strongly expressed.
The model proposed in [10] provides synthetic data fairly close to the characteristics of true data. In this model, the level of expression of a gene is obtained by the superposition of several components.
These components allow one to define the DE genes and the overall level of variability in the data. However, it lacks flexibility in choosing the number of DE genes and the number of genes over-and under-regulated. In this paper, we are interested in generating synthetic microarray data associated with two biological conditions. These conditions correspond to comparison of wild-type samples versus knock-out ones, treated samples versus non-treated ones, etc. For these situations, we use the generic terms of control versus test samples. Two condition biological data can be obtained using either one-color or two-color microarray technology leading to log 2 intensities or log 2 ratios, respectively. We assume that the same reference sample was used for log 2 ratio data.
This paper is organized as follows. In the next section, we present the model used and describe its parameter settings. In Section 3, we present commented results obtained using our model. Conclusions are drawn in Section 4.

Methods
We propose here a model that does not require knowledge of the physical laws governing gene expression to produce data with similar characteristics to the data commonly produced by current platforms. The user provides a subset of parameters that control the behavior of the data generated. We assume the following characteristics for microarray data: 1. There are usually more genes with low intensities than genes with high intensities. One possible explanation for this observation is a non-response of all genes to given biological conditions. Exceptionally different values will be due to technical problems. 4. The total number of DE genes depends on the biological conditions used. The number of over-regulated and under-regulated genes may be different. 5. The variability observed for weakly expressed genes is larger than that observed for highly expressed genes. This is mitigated with today's scanners and tools. However, genes which are not really expressed still have non-zero values, and strongly expressed genes (saturated) receive a single value corresponding to the maximum level of quantification.
In the model described below, the expression levels are log 2 intensities. Ratios are derived from these values.

Model Used
Given n genes, m 1 and m 2 control and test samples, respectively, we used the following model: Values associated with technical problems are cast in vector t i , which has only a few non-zero components (samples) for some genes. The current implementation of our model does not include any technical problem terms, t i , but we intend to take into account this component in the future.
To generate data using model (1) that satisfies the issues given above, we proceed as follows. The beta distribution is used to obtain n values varying between zero and one. Shape parameters (shape 1 and shape 2 ) of this distribution were chosen to produce more small values than higher ones (issue 1).
The n values obtained are scaled to fit real data (issue 2), which vary between a lower bound (lb) and an upper bound (ub). The scaled values represent average expression levels for the n genes (issue 3). Optionally, real data can be used as seed at this step. We assume that intensities for each gene are uniformly distributed around its average level. The range of variation of values from this distribution is assumed to be dependent on the gene average level and is expressed as a percentage through a parameter, α = λ 1 e −λ 1zi , where λ 1 is a parameter setting andz i is the average level of gene, i (issue 5). This leads to vector a i . Two parameters are used for generating vector s i , the percentage of DE genes (pde) and a setting determining the number of up-and down-regulated genes (sym). The first m 1 values (control samples) of s i are set to zero. For gene i, an integer is randomly drawn from the set {1, . . . , 100} and is compared to 100pde to decide whether the corresponding gene is DE. If the answer is yes, a second integer is drawn from the set {1, . . . , 100} and is compared to 100sym to decide up and down assignment. m 2 normally distributed values are generated using a mean (µ de ) and a standard deviation (σ de ). These values (or their opposite) form the second part (test samples) of vector s i for the up-(down-) regulated DE genes. Independent noise vector n i values are obtained using a normal distribution with zero mean and parameter (σ n ) for standard deviation. The details of this algorithm are given in Algorithm 1 and described in the following sections.

Data Size: n
This is the number of probes (genes) on the arrays. It is on the order of thousands, and the default setting is n = 10, 000.

Number of Samples: m 1 and m 2
These parameters are the number of control and test samples, respectively. Typically, the minimum number of samples per condition should be three to allow the use of statistical tests. The maximum value of these parameters rarely exceeds 100 (samples per biological condition); the default setting is seven for both.

Expression Level or Ratio: Ratio
This parameter option generates log 2 intensities or log 2 ratios data. With ratio = 0 (default setting), log 2 intensities data will be generated; otherwise, log 2 ratios data are returned.

Algorithm 1:
Steps of the micorarray data generation model.
1. Generate n values from a beta distribution using shape 1 = 2 and shape 2 as shape parameters, we obtain values z = {z i , i = 1, . . . , n} 2. Transform values of z, such that they vary between settings lb and ub:z = lb + ub × z 3. For each value ofz, generate m 1 + m 2 + 1 uniformly distributed values: and λ 1 is a user setting. Then, the first m 1 + 1 values are used for a vector r 1 , and the last m 2 values are used for a vector r 2 .
• Generate m 2 normally distributed values and add them to r 2 : where λ 2 , µ min de and σ 2 de are settings iii. Else • Generate m 2 normally distributed values and subtract them from r 2 : where λ 2 , µ min de and σ 2 de are settings iv. Set: To obtain more small values than high ones, we use a beta density distribution. Based on many histogram plots using various shape parameter values, we propose to set parameter shape 1 to two and allow the user to choose a value for the parameter shape 2 in the interval [4,8]. The default setting for shape 2 is four.

Log 2 Intensities Variation Range: lb and ub
These parameters specify the log 2 intensities variation range. Quantification of microarray images is usually performed with 16to 20-bit base, leading to 2 16 to 2 20 levels of gray. We suggest using values for parameters lb and ub from intervals [2,6] and [8,16], respectively. These values are typically observed for actual Affymetrix GeneChip c ⃝ array data for gene expression profiling. Default settings are lb = 4 and ub = 14. Observed minimum and maximum expression values will not exactly match the settings because of the variations used and possible additive noise.

Percentage of DE Genes: Pde
This parameter controls the number of differentially expressed genes the user would like to have in the data set. Its values are taken from the interval [0, 1]. A value pde = 0.02 (default setting) means that 2% of the n probes in the data set should be DE.

Number of Up-and Down-Regulated Genes: Sym
This parameter results in nearly the same number of up-and down-regulated probes when sym = 0.5 (default setting). If the value of this parameter is less (greater) than 0.5, we will have more up-(down-) regulated probes in the data set.

Gene Average Level Variation Range: λ 1 (Lambda1)
We assume that the values of each probe are uniformly distributed around an average value. An exponential distribution is used, α = λ 1 e −λ 1zi , where λ 1 is a user setting for decreasing rate. Then parameter, α, controls the width of the uniform distribution and is expressed as a percentage of the average level. λ 1 allows a high variability in weakly expressed genes and, at the same time, a low variability for strongly expressed genes. We use the default value, λ 1 = 0.13, which allows an α ≈ 10% for an average expression level equal to two and an α ≈ 3.54% for an average expression level equal to 10. Increasing λ 1 will lead to more variability for weakly expressed genes and a small variability for strongly expressed genes. Small values for λ 1 (≈0.01) will lead to the same variability for all genes independently of their expression level. For a gene, the fold change is a shift of average expression levels between test and control samples. Assuming that (control and test) sample values of a given gene are uniformly distributed around an average level, the shift comes from a superposition of additional values to the test samples, see Figure 1. We assume a normal distribution for the shift values N (µ de , σ de ). However, the same mean µ de is not used for all DE genes. Hence, we use a minimum value (setting µ min de ) and the exponential distribution to get µ de = µ min de + {λ 2 e −λ 2 }, where λ 2 is another setting. Default settings for these parameters are: µ min de = 1.0, λ 2 = 2 and σ de = 0.5. A higher λ 2 value will lead to a small number of genes having a shift greater than µ min de ; a small λ 2 value leads to the opposite situation. Parameters µ min de and σ de may be chosen using a one sample Student t-test analysis. The statistic of the shift value is z s =   This parameter represents a normal distribution standard deviation for additive noise. The default value is σ n = 0.4. A zero value for this parameter will lead to noise-free data. Too high σ n values can lead to a number of DE genes different to those specified in pde.

Computer Random Generator Seed: Rseed
This parameter is used for computer random number initialization. It will allow one to generate the same data at different times. The default value is 50.

Results and Discussion
To evaluate the performance of the proposed model, we performed simulations, in which we studied the influence of different parameter settings. Their default values are: n = 10, 000, m 1 = m 2 = 7, lb = 4, ub = 14, λ 1 = 0.13, λ 2 = 2, µ min de = 1.0, σ de = 0.5, ratio = 0, sdn = 0.4, shape 2 = 4, pde = 0.02, sym = 0.5). We performed 100 independent simulations by changing the initialization of the generator through parameter, rseed. For each simulation, we performed a Student t-test and selected genes with a p-value less than 0.006. Using the DE information, selected genes were split into two: true and false DE genes. The p-value threshold, 0.006, leads to an expected error (false discovery rate) equal to (0.006 × n)/(pde × n) = 30% for default settings. When studying one parameter, the others are set to their default value.

Parameter Sym
We used the following three values: 0.3, 0.5 and 0.7. For these values, the expected numbers for pairs of down-and up-regulated genes are, respectively, (60, 140), (100, 100) and (140, 60). Figure 3 shows the boxplot of results obtained. The median numbers of true down-and up-regulated genes observed are (46, 107), (75.5, 76.5) and (107, 46) for the above values of parameter sym, respectively. Hence, the recovery powers of the Student t-test are 76.5%, 76% and 76.5%.  The median number of false DE genes is 49 for the three values of parameter sym.

Parameter σ n
We used three values: 0.2, 0.4 and 0.6. Figure 4 shows the boxplots of the results obtained. The median numbers of down-and up-regulated genes observed are (89, 89), (75.5, 76.5) and (57, 58), leading to detection powers of 89%, 76% and 57.5%, respectively. The t-test detection power decreases when σ n increases.   We performed simulations to examine the influence of these parameters. For each parameter, 100 simulations were used, and the results obtained are summarized in Table 1. Increasing the parameter µ min de setting introduces more change for the DE genes, while its decrease leads to the opposite effect. Parameter σ de acts as noise. The effect of the modification of some parameters is investigated further in the MA plot representations described in the following paragraph.

Volcano and MA Plots
Using default settings, we performed one simulation. Then we computed the Student t-test p-value and fold change for all genes. These values were used in the volcano plot of Figure 5. Red circles represent genes having a p-value less than 0.01 and a fold change greater than two or less than 0.5. Intensity measurements of two samples can be used to create two new variables: M = I x 2 − I x 1 and A = 0.5(I x 2 + I x 1 ), where I x 2 and I x 1 are log 2 intensities of samples x 2 and x 1 , respectively. A value, one (−1) for M means that the corresponding gene is up-(down-) regulated two-fold. A plot of M (log 2 ratio) versus A (log 2 intensities average) is denoted "MA plot" [11]. The MA plots in Figure 6 are obtained using either two control samples (panel A) or one control and one test sample (panel B).
Additional MA plots in Figures 7-10 showing the effect of some parameter settings. A small value for λ 1 leads to a less dense cloud of points. A larger change is observed for higher values of λ 2 than for smaller ones. The same applies to parameter, σ de . Increasing the parameter shape 2 value leads to a decrease of the dynamic range of the data.

Discussion
The choice of some setting parameters for the proposed model is easy and can be dictated by the experimental design. This applies to n, m 1 , m 2 and ratio. Parameters lb and ub and sym and pde concern the dynamic range of variation of the data to generate and define the DE genes. The intervals indicated for lb and ub are those observed for data from common platforms. Parameters, shape 2 , lb and ub, have no effect if real microarray data are used as a seed. The number of DE genes (pde) and the proportion of under-and over-regulated genes (sym) is at the discretion of the user. The parameter, rseed, allows one to produce the same data at different times using the same computer. This parameter is also useful for generating test data for different analysis algorithms. Settings, λ 1 , λ 2 , µ min de , σ de and σ n , control the global behavior of the data generated. More precisely, λ 1 allows one to make gene changes dependent on the average level of expression. λ 2 introduces variation in the expression changes for DE genes. These changes are defined by µ min de and σ de . Parameter, σ de , also acts as noise. The parameter, σ n , allows one to perturb the data generated. The example of Figure 4 shows that the data obtained can differ from those expected for large values of σ n . The signal to noise ratio (defined as the ratio of standard deviations) for default settings is 6.25. This rises to 4.16 for σ n = 0.6.
An interesting microarray study was reported in [12]. In that study, the same RNA samples were processed by many laboratories using three leading microarray platforms: Affymetrix (five labs), two-color cDNA (three labs) and two-color oligonucleotides (two labs). The results presented show a good agreement across-platforms in contrast to some results previously reported in the literature, see, for instance, references cited in [12]. The microarray data generated for the study described in [12] are available from the Gene Expression Omnibus [13] under accession number GSE2521. These data can be used as seed for our model, which can then be integrated into user friendly data analysis software, such as Partek Genomics Suite, GeneSpring GX, etc., for demonstration and/or teaching purposes.

R Code
For immediate use of the proposed model, we provide an R code function madsim.R (MicroArray Data Simulation Model), which is deposed as a package on the Comprehensible R Archive Network (CRAN) server for download [14]. The outputs of this function are the data generated (xdata) and the indexes of DE genes (xid). Real data can be used as seed for each gene. An example of such data is available in the data folder of the package. Further explanations are available from the package's help function.

Conclusions
We proposed in this paper a simulation model of microarray data. This model is very flexible and allows one to generate data with similar characteristics to the data commonly produced by current platforms. We showed a commented example of its possible use. We considered the case of data from two biological conditions. This model can be extended to multiple biological conditions in different ways: (a) modify to take into account additional biological conditions and several levels for the parameters, µ de and σ de , and (b) use as is, then place data side by side.