Agn Ès Hamon

Marie Tosolini

Bernard Ycart
email: bernard.ycart@imag.fr

Fr Éd Éric Pont

Jean-Jacques Fourni

Simultaneous growth of two cancer cell lines Hamon et al Simultaneous growth of two cancer cell lines evidences variability in growth rates

Keywords:

Cancer cells co-cultured in vitro reveal unexpected differential growth rates that classical exponential growth models cannot account for. Two noninteracting cell lines were grown in the same culture, and counts of each species were recorded at periodic times. The relative growth of population ratios was found to depend on the initial proportion, in contradiction with the traditional exponential growth model. The proposed explanation is the variability of growth rates for clones inside the same cell line. This leads to a log-quadratic growth model that provides both a theoretical explanation to the phenomenon that was observed, and a better fit to our growth data.

Introduction

Since emergence of resistant cells underlies time to relapse for cancer patients undergoing chemotherapy, the growth rate of these tumor cells is a crucial issue. Cancer therapies usually yield undetectable levels of residual and resistant cancer stem cells (CSC) in patients. Upon repeated mitosis however, CSC can seed a cell progeny that progressively reconstitutes tumors, but the proportion and mitotic rate of such CSC are highly variable in treated patients. The classical exponential growth model predicts that the relative growth of the fast-growing clones should increase exponentially with time, regardless of their initial rates in patients. On the other hand however, this model is challenged by heterogeneity of the clonal progeny from a cancer cell and the resulting Darwinian selection in this progeny for access to nutriments. [START_REF] Chevin | On measuring selection in experimental evolution[END_REF][START_REF] Manna | The high-throughput yeast deletion fitness data and the theories of dominance[END_REF] To investigate this, we grew either separately or together, two different non-interacting human cancer cell lines, in cell cultures with unlimited medium supply, and modeled the cell growth rates observed in the co-cultures. The exponential growth model is so elementary and has been known for such a long time, [START_REF] Powell | Growth rate and generation time of bacteria with special reference to continuous culture[END_REF][START_REF] Kimmel | Branching processes in Biology[END_REF] that it seems almost too simple to actually fit real cell growth data. [START_REF] Koutsoumanis | Stochasticity in colonial growth dynamics of individual bacterial cells[END_REF] Yet, for a given cell line grown in unlimited supporting medium, an excellent linear fit is usually observed for the logarithm of population size against time. [START_REF] Deenick | Stochastic model of T-cell proliferation: a calculus revealing IL-2 regulation of precursor frequencies, cell cycle time, and survival[END_REF][START_REF] Zilman | Stochastic models of lymphocyte proliferation and death[END_REF] Our experiments were conducted with two well known lab strains: RL (non-Hodgkin's lymphoma B cell line: ATCC CRL-2261 [START_REF] Beckwith | Phorbol ester-induced, cell-cycle-specific, growth inhibition of human blymphoma cell lines[END_REF] ) and THP-1 (cell line derived from an acute monocytic leukemia patient: TCC TIB-202 [START_REF] Tsuchiya | Establishment and characterization of a human acute monocytic leukemia cell line (THP-1)[END_REF] ). As a control experiment, the two cell lines were grown separately. An excellent log-linear fit was observed. Then, both strains were grown in the same solution. No interaction between the two species could occur, other than possible competition for nutrient from culture medium. This was avoided again by maintaining a sufficient supply of medium by volume, and continuous renewal. The initial proportions of the faster growing strain RL were fixed at 0.5%, 1%, 5%, and three replicates were made for each initial proportion. With the classical exponential growth model, the relative proportion of RL vs. THP-1 would be predicted to increase exponentially in time, at a rate which is independent from the initial proportion. Somewhat unexpectedly, this turned out to be false. Figure 1 presents a plot of the ratio of observed RL vs THP-1 counts on a logarithmic scale. The time scales have been shifted so that the origin corresponds to the time at which each proportion reaches 5%. The slope of the regression line decreases as the initial proportion of the faster strain increases: the slope with an initial proportion 5% (red) is smaller than that with initial proportions 1% (green) and 0.5% (blue).

In vitro experiments with simulateous growth of two or more microorganisms have long been carried through: see Dykhuizen [START_REF] Dykhuizen | Experimental studies of natural selection in bacteria[END_REF] for a review. The variability of growth rates in human leukemia cell clones has been studied by Tomelleri et al. [START_REF] Tomelleri | A quantitative study of growth variability of tumour cell clones in vitro[END_REF] However, to the best of our knowledge, this is the first instance describing an experiment with two different cancer lines, and showing the phenomenon of dependence of the ratio growth rate on the initial condition, illustrated on Figure 1.

The objective of this paper is to propose a stochastic growth model explaining the phenomenon, and show that the fit of the data by that model is better than with the exponential growth model.

It has long been known that exponential proliferation is a valid approximation, only on a certain fraction of the observation time. [START_REF] Novick | Growth of bacteria[END_REF] Many different models have been proposed as growth curves. [START_REF] Zwietering | Modeling of the bacterial growth curve[END_REF][START_REF] Peleg | Microbial growth curves: what the models tell us and what they cannot[END_REF] At the beginning of a cell growth experiment, a lag phase [START_REF] Brook | Variability in the cell cycle and control of proliferation[END_REF][START_REF] Delignette-Muller | Relation between the generation time and the lag time of bacterial growth kinetics[END_REF][START_REF] Baty | Estimating the bacterial lag time: which model, which precision[END_REF] is usually observed; this is the case in our data.

The lag phase could partly account for the phenomenon investigated here. Indeed, when starting from a 0.5% proportion, the lag phase has elapsed when reaching 1%, but if one already starts with a 1% proportion, the lag phase only begins. Yet the lag phase does not explain differential growth rates after all cultures have reached a 5% proportion, since at that time, the lag phase has elapsed in the first two cases. Another simple explanation is proposed here: intrinsic variability of growth rates. [START_REF] Koutsoumanis | Stochasticity in colonial growth dynamics of individual bacterial cells[END_REF][START_REF] Tomelleri | A quantitative study of growth variability of tumour cell clones in vitro[END_REF] Here, the notion of growth rate is understood in the sense of branching processes, [START_REF] Kimmel | Branching processes in Biology[END_REF] as a "large scale approximation" that applies to the whole clone stemming from a given cell, and not just to The three replicates for each initial condition are marked by circles, triangles and squares. All curves have been shifted to the first day where the proportion RL/THP-1 passes 5%. Red marks correspond to an initial proportion of 5%, green marks to 1%, blue marks to 0.5%. The three regression lines are represented with corresponding colors. The slopes (time unit: hour) are 0.0084 (red line, initial proportion 5%), 0.0138 (green line, initial proportion 1%), and 0.0154 (blue line, initial proportion 0.5%). that cell. Therefore, the variation of growth rates can only be genetic: to each cell present at the beginning of the experiment is associated one value, which will be the growth rate of the whole clone stemming from that cell. If growth rates among RL vary, the proportion of fast dividing mitotic cells among RL cell clones will gradually increase. When reaching the proportion of 5%, there will be more fast breeders among RL if the initial proportion was 0.5%, than if it was 1%. This intuitively explains why the estimated growth rate over a given time interval, is larger when the proportion of RL reaches 5%, after starting from 0.5%.

Mathematically, it will be shown that assuming variable growth rates among the cells of a same species, naturally leads to a log-quadratic model on the population growth, instead of the traditional log-linear (or exponential) model. It will be shown that the log-quadratic model induces a better fit on our data.

Using that model, the observed phenomenon can be explained and quantified. Indeed, if the actual growth is log-quadratic instead of log-linear, the fit of a log-linear model yields estimated slopes that vary with the initial condition: Figure 2 illustrates that theoretical explanation. The derivation of the log-quadratic model from the hypothesis of variable growth rates uses the cumulant generating function of the random growth rates. A similar explanation had been proposed by Hansen. 18

Material and methods

Experimental methods

The cell line THP-1 (ATCC TIB-202) derives from human acute monocytic leukemia. It has a monocyte morphology and expresses the cell surface marker CD13. The cell line RL (ATCC CRL-2261) was derived from human non-Hodgkin's lymphoma and expresses the cell surface marker CD20. These two cell lines were cultured as indicated by the supplier (ATCC www.lgcstandardsatcc.org) at 37 o C and 5% CO2 in liquid medium RPMI-1640 (LONZA, Levallois, France) supplemented with 10% heat-inactivated fetal calf serum (FCS), 2mM L-glutamine, 100 U/ml penicillin and 100g/ml streptomycin (Invitrogen, Cergy Pontoise, France). This medium contains inorganic salt, amino acid, vitamins and D-glucose (2g/L). THP-1 and RL, alone and in competition, were cultured in T75 flasks with 50mL of medium without agitation. At the beginning of the culture, for all the conditions, cells were seeded at 0. the differents cell cultures were counted and if necessary, diluted with complete medium if the cell concentration was higher than 0.7 × 10 6 cells/mL to adjust at 0.3 × 10 6 cells/mL. These concentrations lead to no competition for nutrient. The dilutions of the culture were adjusted to the cell growth in each flask: RL cell culture was more diluted than THP-1 cell culture for example.

Cells from competition cultures were also analysed by flow cytometry to determine the percentage of each cell line in the culture: Cells were centrifuged, washed with PBS and incubated 10 min with antibodies against CD20 coupled with the fluorochrome APC-Cy7 to identify RL cells, and against CD13 coupled with the fluorochrome PE (BD Biosciences, Pont de Claix, France) to identify THP-1 cells. Cells from single cell line culture were used as controls. The fluorescence of 50000 cells was then analysed using a BD LSR II cytometer. Two sets of experiments were carried through 35 days, with daily measurements. In the first set, the two cell lines RL and THP-1 were grown in separate culture flasks, in duplicate. In the second set, the two cell lines were grown in the same culture flask, in triplicate, with unlimited amount of nutrient in each case. Three initial proportions of RL (the faster growing strain) were considered: 0.5%, 1%, and 5%. For each set of experiment, each day of culture, and each replicate, numbers of cells of each type were recorded. The dataset is available upon request. On these data, different least square fits of the logquadratic model (6) were performed, for each of the two separate growths (first set of experiments), and for simultaneous growth (second set).

Mathematical model

In this section, a mathematical derivation of a log-quadratic growth model, based on variable growth rates, is proposed.

Consider first the classical model of exponential growth for a single clone, stemming from the general theory of branching processes. [START_REF] Harris | The theory of branching processes[END_REF][START_REF] Athreya | Branching processes[END_REF][START_REF] Kimmel | Branching processes in Biology[END_REF] From a single cell at time 0, the clone grows to size N (t) at time t. Under fairly general hypotheses on the division time distribution, there exists a positive constant b, the growth rate (also called Malthusian parameter), such that almost surely:

lim t→+∞ e -bt N (t) = C , (1)
where C is a random variable with finite expectation and variance. This is one of the basic results of the theory of branching processes. [START_REF] Harris | The theory of branching processes[END_REF][START_REF] Athreya | Branching processes[END_REF] Thus it is reasonable to assume N (t) = Ce bt as a model of growth curve for a single clone. Assume now that the population grows from a large number n of identical initial cells. For i = 1, . . . , n, let N i (t) be the size at time t of the clone stemming from cell i:

N i (t) = e bt C i , (2)
where the C i are independent identically distributed random variables. The total population at time t is:

N (t) = n i=1 N i (t) .
By the law of large numbers, almost surely:

lim n→∞ N (t) n = e νt E(C) ,
where E(C) denotes the mathematical expectation of the random variable C.

This justifies the classical log-linear model:

log(N (t)) = a + bt , (3)
where a = log(N (0)). General references on log-linear models are Mair [START_REF] Mair | Interpreting standard and non standard log-linear models[END_REF] and von Eye and Mun. [START_REF] Eye | Log-linear modeling: concepts interpretation, and application[END_REF] Consider now a second population growing according to the same model, and denote by M (t) its size at time t.

log(M (t)) = a ′ + b ′ t .
Assume b > b ′ (the first population grows faster). Then the ratio R(t) = N (t)/M (t) also follows a log-linear model.

log (R(t)) = (a -a ′ ) + (b -b ′ )t .
Whatever the interval of time it is observed on, the growth rate bb ′ does not depend on the interval nor on the initial proportion. This is contradicted by our observation (Figure 1). Assume now that clones stemming from different initial cells may have different growth rates. The new model is:

N i (t) = e Bit C i , (4)
where (B i , C i ) are independent and identically distributed copies of a random couple (B, C). The joint distribution of (B, C) is of course unknown, and we shall make the technical assumptions that B and C are independent, and that B has faster than exponential decaying tails. By the same argument of law of large numbers, the global population N (t) should satisfy:

log(N (t)) = a + log E e Bt .
Note that the function log E e Bt exists for all t 0 if the distribution of B has faster than exponential decaying tails. This function is the cumulant generating function of B, [START_REF] Kendall | The advanced theory of statistics: distribution theory[END_REF] well known in large deviation theory. [START_REF] Dembo | Large deviations techniques and applications[END_REF] Let µ be the expectation of B, σ its standard deviation and γ 1 its skewness. Then the first three terms of the Taylor expansion of log E e Bt are:

log E[e Bt ] = µ t + σ 2 2 t 2 + σ 3 γ 1 6 t 3 + o(t 3 ) , (5)
In the particular case where B follows the Gaussian distribution, the first two terms give the exact expression:

log(E[e Bt ]) = µ t + σ 2 2 t 2 .
In that case, the growth of N (t) is quadratic in logarithmic scale:

log(N (t)) = a + bt + ct 2 , (6)
with a = E[N 0 ], b = µ, and c = σ 2 /2. Equation ( 6) will be referred to as logquadratic model: see Chapter 9 of von Eye and Mun, [START_REF] Eye | Log-linear modeling: concepts interpretation, and application[END_REF] and Stone et al 25 for an application in a similar context. Assuming that the distribution of B is Gaussian may seem unrealistic, but whatever the distribution of B, if its expectation is µ = b and variance σ 2 = 2c, equation ( 6) remains true as a second order approximation, because of (5). This justifies the use of ( 6) as a model, in case of variable growth rates.

If two populations grow according to a log-quadratic model, then the ratio of the two population sizes does too. Denote again that ratio by R(t), assuming that the choice has been made to put the faster growing population on the numerator, so that R(t) increases.

log(R(t)) = a + bt + ct 2 , (7)
with a = log(R(0)) and b, c > 0. In practice, growth rates are estimated by a log-linear regression over a given interval, say [T 1 , T 2 ]. This amounts to approximating (7) by: log(R(t)) = α + βt , where α and β are optimal in the sense of mean squares:

(α, β) = arg min T2 T1 (a + bt + ct 2 -α -βt) 2 dt . (8)
The solution of ( 8) is easily obtained:

α = a - c 6 (T 2 1 + T 2 2 + 4T 1 T 2 ) and β = b + c(T 1 + T 2 ) .
For a fixed span T 2 -T 1 , the "equivalent growth rate" β increases as T 1 increases (see Figure 2 for an illustration). This explains the phenomenon evidenced by Figure 1. More precisely, let T 1 be the time at which R reaches the value R 1 > R(0):

T 1 = 1 2 -b + b 2 + 4c log R 1 R(0) .
The equivalent growth rate on a time interval of duration t after T 1 will be β = b + c(2T 1 + t). It will be larger than the growth rate on an interval of same width starting at 0, which is b + ct.

Thus the log-quadratic model ( 6) provides a theoretical explanation for the phenomenon of differential growth rates, that has been observed. As will be shown in the next section, it also provides a better fit to our data.

Results

Separate growth of RL

Let Y ik denote the logarithm of cell count at time t k and replicate i (i = 1, 2).

We consider the following model :

Y ik = a + bt k + ct 2 k + ε ik , (9)
where ε ik are centered Gaussian random variables with common standarddeviation. For RL cells, it turned out that the coefficient c, that we shall call "curvature", was not significantly different from zero (P = 0.698). Therefore a linear model without quadratic term was fitted. fitted model is excellent (R 2 = 0.99). The QQ-plot of residuals (figure 3) is close to linear, the plot of residuals vs. time does not show any mis-specification of the non-random part nor heteroscedasticity problem. The Durbin-Watson test (P = 0.12) and the runs test (P = 0.69) indicate no violation of the hypothesis of error independence.

Separate growth of THP-1

The model remains the same, see equation (9). At first, it was fitted to the full data set. Three observations at the end of the experiment period, were detected as outliers, and therefore excluded from the final analysis. Table 2 reports the estimated coefficients. As expected, the growth rate of THP-1 is significantly smaller than that of RL (P < 0.0001).

ACF of residuals

ACF of residuals

Simultaneous growth of RL and THP-1

Let Y ijk denote the logarithm of the ratio of RL by THP-1 cell counts at time t k , where i denotes the replicate (i = 1, 2, 3) and j the initial nominal value of the ratio, that will be called dilution. Indices j = 1, 2, 3 correspond to dilutions 0.5%, 1%, 5%. Four different models were considered:

Y ijk = a + a j + bt k + ǫ ijk (M0) Y ijk = a + a j + bt k + b j t k + ǫ ijk (M1) Y ijk = a + a j + bt k + b j t k + ct 2 k + ǫ ijk (M2b)
Y ijk = a + a j + bt k + ct 2 k + c j t 2 k + ǫ ijk (M2c)
Y ijk = a + a j + bt k + b j t k + ct 2 k + c j t 2 k + ǫ ijk (M3)
Model (M0) is the simplest: the expected log ratio Y ijk is modeled by a straight line, the slope b of which does not depend on dilution. In model (M1), the slopes b + b j may depend on dilution. In the actual fit, the slopes are found to decrease as the initial proportion of RL increases. This is coherent with Figure 1. However, it hides the relevance of the quadratic models. Model (M3) is the complete log-quadratic model: both the slopes b + b j and the curvatures c + c j may depend on dilution. Models (M2b) and (M2c) are embedded into (M3): in model (M2b) the slope does not depend on dilution and in (M2c) the curvature does not depend on dilution. For all four models, the linear fit was computed, then pairs of embedded models were tested by the Fisher test of analysis of variance. The results are presented in Table 3: the degrees of freedom df, the Fisher test statistic F, and the significance p-values are given. The conclusions are the following. The first three comparisons are significant, ie the bigger model is better than the embedded one. The (M3) vs. (M2c) comparison is not. The conclusion of the four comparisons is that the best fitted model is (M2c). This indicates that if curvatures are included in the model, the slopes do not significantly depend on dilution. This is coherent with the theoretical derivation of the log-quadratic model (7). In model (M2c) the estimated slope is b = 6.1 × 10 -3 and 95% confidence interval on b is [5.6 × 10 -3 ; 6.6 × 10 -3 ]. Recall from ( 7) that the slope b of (M2c) should be understood as the difference between the slopes of models ( 9) for RL and THP-1. From the two previous sections, the estimated difference is 1.01 × 10 -2 , which is above the confidence interval on b in (M2c).

Next, we tested the three pairwise differences of curvatures c j in the accepted model (M2c). The results are presented in Student test statistic t, and the p-value are given (degrees of freedom: 293). All three differences are significant. The estimated values of the curvatures c + c j 's are given in Table 5. It turns out that c 1 > c 2 > c 3 . So the curvature c + c j decreases as the initial proportion of RL increases. This is a similar phenomenon as was observed on Figure (1). Indeed, when curvatures are neglected (model (M1)), the slopes were found to be decreasing as the initial proportion of RL increases. The theoretical explanation is given by the mathematical model. Model (M1) amounts to keeping only the first term in the Taylor expansion (5). Since the next term is positive, the adjusted values of the slopes increase with time. Model (M2c) considers the first two terms in (5), neglecting the third one. If that neglected term is positive, then the same effect will occur: ajusted values of curvatures increase with time. The third term is proportional to the skewness of B. Observing decreasing curvatures as in Table 5 is an indication that the skewness of B may be positive, where B is the (random) difference in growth rate between RL and THP-1. Koutsoumanis and Lianou [START_REF] Koutsoumanis | Stochasticity in colonial growth dynamics of individual bacterial cells[END_REF] have proposed a logistic distribution as a model for growth rate variability. That distribution has null skewness. We conjecture that distributions with positive skewness provide better models for variable growth rates.

Conclusion

That unchecked populations grow exponentially fast is a well known fact, backed up by countless experiments, which validate the mathematical theory Dilution j c + c j j = 1 (0.5%) 5.3 × 10 -6 j = 2 (1%)

4.96 × 10 -6 j = 3 (5%)

1.6 × 10 -6 of branching processes. [START_REF] Kimmel | Branching processes in Biology[END_REF] To any clone stemming from a single cell can be associated an exponential growth rate, also called Malthusian parameter. It can be seen as the slope over a large interval of time, of the line fitting logarithms of the number of cells against time. What is questioned here, is the idea that clones stemming from different cells in a given strain, should have the same growth rate. Unlike Koutsoumanis and Lianou 5 or Tomelleri et al, [START_REF] Tomelleri | A quantitative study of growth variability of tumour cell clones in vitro[END_REF] we do not provide direct evidence for the intrinsic variablility of growth rates, but instead an indirect proof, coming from a simultaneous growth experiment. [START_REF] Dykhuizen | Experimental studies of natural selection in bacteria[END_REF] It consisted in growing in the same vessels two cancer cell lines, RL [START_REF] Beckwith | Phorbol ester-induced, cell-cycle-specific, growth inhibition of human blymphoma cell lines[END_REF] and THP-1. [START_REF] Tsuchiya | Establishment and characterization of a human acute monocytic leukemia cell line (THP-1)[END_REF] If there existed a single growth rate for all RL clones, and another for all THP-1 clones, then the ratio should grow exponentially, the rate being the difference of the two growth rates. In that case, the growth rate of the ratio should not depend on the initial proportion of RL vs. THP-1. Our observations disproved this: the growth rate of the ratio was found to increase as the initial proportion of RL decreased (Figure 1). Assuming that growth rates may vary among clones provides both intuitive, and theoretical explanations.

The intuitive explanation is the following. Consider a growth rate as attached to each cell of a given clone. If clones may grow at different rates, the proportion of cells in faster growing clones will gradually increase. In other words, the distribution of growth rates at increasing times will be shifted toward larger values. This explains why, when the initial proportion of RL cells is 0.5%, at the time it reaches 5%, the population of RL contains more fast breeders than at time 0. Therefore, the (apparent) growth rate for an initial proportion of 0.5% is larger than for an initial proportion of 5%. The theoretical explanation is the following. If growth rates of different clones are considered as independent random variables with a positive variance, then the model fitting the logarithms of cell numbers against time must contain a quadratic term, proportional to the variance of growth rates: variable growth rates imply that higher order terms must be added to the classical loglinear model. Now if a population grows according to a log-quadratic model, and a log-linear model is fitted instead, then the estimated slope over an interval of time should increase as the interval moves to the right (see Figure 2). This will overestimate the mean growth rate.

To validate our theoretical explanation, we had to compare the fits of the log-linear and log-quadratic models on our experimental data. On separate growth data, the log-linear model was better than the log-quadratic model on RL, the contrary was true for THP-1. On simultaneous growth data, the logquadratic clearly provided a better fit. This may seem paradoxical: indeed, the same model should be adopted for separate and simultaneous growths. The explanation of this apparent contradiction is statistical. The estimated curvature terms are in all cases smaller by several orders of magnitude than the estimated slopes. Therefore, the log-linear and log-quadratic models can hardly be distinguished when the cell counts range over several orders of magnitude, as in separate growths. This cannot be the case on simultaneous growth data, where the ratios range from a few percents to 100%. We believe that, on a separate growth experiment, if more values were collected at the beginning, then the log-quadratic model would provide a better fit.

There remains the issue of a probabilistic model to be fit on variable growth rates. Our derivation of the cumulant generating function shows that classical models of positive random variables, such as Gamma, Log-normal, or Logistic distributions 5 cannot be used here. Indeed their exponentially decaying tail implies that the equivalent growth rate would become infinite at finite time, which is not realistic. So a truncated model would have to be used instead. In any case, there would remain to adjust the chosen distribution to actual data. Ideally, these data should be collected from the observation of colonial growth of individual cells, such as reported by Koutsoumanis and Lianou for Salmonella enterica, [START_REF] Koutsoumanis | Stochasticity in colonial growth dynamics of individual bacterial cells[END_REF] or Tomelleri et al 11 on leukemia cells. This will be the object of future work.

Figure 1 :

 1 Figure1: Logarithms of ratios of RL/THP-1. The three replicates for each initial condition are marked by circles, triangles and squares. All curves have been shifted to the first day where the proportion RL/THP-1 passes 5%. Red marks correspond to an initial proportion of 5%, green marks to 1%, blue marks to 0.5%. The three regression lines are represented with corresponding colors. The slopes (time unit: hour) are 0.0084 (red line, initial proportion 5%), 0.0138 (green line, initial proportion 1%), and 0.0154 (blue line, initial proportion 0.5%).

3 ×Figure 2 :

 32 Figure 2: Theoretical explanation for observed differential growth rates. The figure represents a quadratic growth in time, at 3 different intercepts: log(0.05) (red), log(0.01) (green), log(0.005) (blue). The three solid curves are parabolas. The green and blue dashed lines are linear fits over an interval starting at the point where the corresponding parabola reaches log(0.05). The green and blue solid lines are time shifts, illustrating the differential slopes.

Figure 4

 4 presents the residual analysis. Contrarily to the RL case, the curvature c is significantly positive (P = 8.6 × 10 -11 ). The proportion of the variation of Y ik explained by the fitted model is excellent (R 2 = 0.99). The 95% confidence interval of the mean THP-1 growth rate b is [0.0204; 0.0227]. This corresponds to a doubling time between 30.5 and 34 hours, i.e. slightly below the values given in Tsuchiya et al[START_REF] Tsuchiya | Establishment and characterization of a human acute monocytic leukemia cell line (THP-1)[END_REF] (35 to 50 hours), and above those of Tsuchiya et al 26 (24 to 30 hours).The QQ-plot of residuals (figure4) is close to linear, the plot of residuals vs. time does not show any mis-specification of the non-random part nor heteroscedasticity problem. The Durbin-Watson test (P = 0.16) and the runs test (P = 0.75) indicate no violation of the hypothesis of error independence.

+

  + + + + + + + + + + + + + + ++ + +

Figure 3 :+

 3 Figure 3: Fit and validation of the log-linear model: RL cells.

Figure 4 :

 4 Figure 4: Fit and validation of the log-quadratic model: THP-1 cells.

Table 1

 1

	reports the estimated

Table 1 :

 1 Estimations for the logarithm of RL cell counts.

	Coefficient Estimate Std. Error	P
	a	12.37	5 × 10 -2	< 2 × 10 -16
	b	0.0212	4 × 10 -4	< 2 × 10 -16
	c 5 × 10 -6	7 × 10 -7	8.6 × 10 -11

Table 2 :

 2 Estimations for the logarithm of THP-1 cell counts.

Table 4 :

 4 the value of the

	Embedded models	df	F	p-value
	(M3) vs. (M2b)	(2, 291) 24.67 1.3 × 10 -10
	(M2b) vs. (M1)	(1, 293) 22.8	8 × 10 -39
	(M1) vs. (M0)	(2, 294)	118	1.8 × 10 -38
	(M3) vs. (M2c)	(2, 291) 2.68	0.0705

Table 3 :

 3 Tests of embedded models for the log ratio of RL vs THP-1 cell counts.

	null hypothesis	t	p-value
	c 1 = c 2	2.4	0.016
	c 1 = c 3	17.9 6.3 × 10 -49
	c 2 = c 3	16.1 4.5 × 10 -42

Table 4 :

 4 Pairwise tests for differences in curvatures c j in model (M2c).

Table 5 :

 5 Estimated curvatures c + c j in model (M2c).

S. Tsuchiya, N. Minegishi, M. Minegishi, T. Sato, and T. Konno. Adaptation of a human monocytic leukemia cell line (THP-1) in a protein-free chemically defined medium. Tohoku. J. Exp. Med, 150(4):455-465, 1986.

Acknowledgements

This work was supported by Laboratoire d'Excellence TOUCAN (Toulouse cancer). The authors are indebted to the anonymous referees for important suggestions.

Disclosure

The authors report no conflicts of interest in this work.


Bernard Ycart
email: bernard.ycart@imag.fr

Konstantina Charmpi
email: konstantina.charmpi@imag.fr

Sophie Rousseaux
email: sophie.rousseaux@ujf-grenoble.fr

Jean-Jacques Fournié

Large scale statistical analysis of GEO datasets

The problem addressed here is that of simultaneous treatment of several gene expression datasets, possibly collected under different experimental conditions and/or platforms. Using robust statistics, a large scale statistical analysis has been conducted over 20 datasets downloaded from the Gene Expression Omnibus repository. The differences between datasets are compared to the variability inside a given dataset. Evidence that meaningful biological information can be extracted by merging different sources is provided.

Background

Many genomewide expression datasets have been published during the past ten years. Repositories, such as the Gene Expression Omnibus (GEO) database [START_REF] Edgar | Gene Expression Omnibus: NCBI gene expression and hybridization array data repository[END_REF], have made available an impressive wealth of data. Using them as a whole, instead of restricting statistical studies to one particular dataset, is tantalizing. Two recently published R/Bioconductor packages [START_REF] Heider | virtualArray: a R/bioconductor package to merge raw data from different microarray platforms[END_REF][START_REF] Taminau | inSilicoMerging: Collection of Merging Techniques for Gene Expression Data[END_REF] provide various tools for merging datasets coming from different studies. However, a serious doubt has been cast by Haibe-Kains et al. [START_REF] Haibe-Kains | Inconsistency in large pharmacogenomic studies[END_REF], after comparing two large scale pharmacogenomic studies: whereas both studies had a good overall correlation, important discordances could be observed. Thus, the following crucial question remains to be answered: is it statistically legitimate to merge datasets coming from different studies? An attempt at answering this question is reported here.

Merging different datasets, requires prior checking that the information they contain is compatible, and hence that detected differences between gene expressions under different conditions are not artifacts, due to experimental or data processing methods. An obvious obstacle to simultaneous treatment is that expression data collected under different experimental conditions and/or platforms usually have incompatible 1 distributions, which differ sometimes by several orders of magnitude [START_REF] Kuo | Analysis of matched mRNA measurements from two different microarray technologies[END_REF][START_REF] Mah | A comparison of oligonucleotide and cDNA-based microarray systems[END_REF]. A solution is provided by robust (or distribution-free) statistics [START_REF] Gibbons | Nonparametric statistical inference[END_REF][START_REF] Héritier | Robust methods in biostatistics[END_REF]. Robust methods amount to replacing actual values by ranks, or equivalently by empirical distribution functions or van der Waerden's normal scores [7, p. 309]. This idea has already been applied to expression data in several papers, including [START_REF] Tsodikov | Adjustments and measures of differential expression for microarray data[END_REF][START_REF] Warnat | Cross-platform analysis of cancer microarray data improves gene expression based classification of phenotypes[END_REF][START_REF] Breitling | Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments[END_REF]. However, to the best of our knowledge, a large scale analysis assessing the reproducibility of information from one dataset to another, is still missing.

We have conducted such an analysis over 20 GEO datasets, totalling 17 745 genomewide expression samples.

For the data treatments presented here, the statistical language R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] has been used. Our set of functions, together with a manual, has been made available online as supplementary information. Throughout the article, we consider data matrices (also called assay data in [START_REF] Gentleman | Bioconductor: open software development for computational biology and bioinformatics[END_REF]) as containing expression data relative to a set of genes. Each row corresponds to a different gene symbol, or feature, each column to a different data vector or sample (see Table 2 in [START_REF] Edgar | Gene Expression Omnibus: NCBI gene expression and hybridization array data repository[END_REF]). Such a matrix is deduced from raw datasets, available on the GEO repository, though standard treatments: annotation and reduction [START_REF] Gentleman | genefilter: genefilter: methods for filtering genes from microarray experiments[END_REF]. Several R packages [START_REF] Davis | GEOquery: a bridge between the Gene Expression Omnibus (GEO) and Bio-Conductor[END_REF][START_REF] Taminau | inSilicoDb: Access to the InSilico Database[END_REF] that perform these operations and output data matrices such as considered here, are available. We have encoded our own functions. We have chosen a data structure in which each data matrix is paired with its information matrix. The columns of the information matrix are labelled by the same numbers as the paired data matrix. Its rows contain the different information fields of the data. Our focus here is on overexpression or underexpression of genes, in different tissues or cancer types.

Our objective was twofold. On the one hand, we wanted to check whether the information on genes, contained in different data matrices, was compatible, and to which extent. This was done on a set of 20 different matrices. Various statistical treatments were performed. The first one consisted in computing correlations between median columns of the matrices. Vectors of pairwise correlations between rows were also compared. Then multivariate analysis over assays of gene symbols was applied: Wilcoxon and Kruskal-Wallis tests, factor and principal component analysis (PCA). The results were compared to those obtained by sorting a single matrix according to different keywords. All comparisons showed not perfect, but highly significant correlations. However, it was also found that in all cases, a sizeable proportion of symbols were good discriminators of the different matrices. But this was also found to hold between two submatrices inside a given dataset. Therefore, it cannot be regarded as an obstacle to merging different datasets. On the other hand, we wanted to know whether biological information could be consistently retrieved from matrices collated from different sources. Two merged sets of matrices were made. The first one came from general cancer cell datasets, from which samples of breast and lung tumors were extracted. The second one was made of blood RNA datasets, coming from healthy individuals, or from leukemias. In both cases, evidence that already known biological information could be extracted from merged matrices was found.

Methods

The datasets available from the GEO repository [START_REF] Edgar | Gene Expression Omnibus: NCBI gene expression and hybridization array data repository[END_REF], collate sets of expression vectors, or samples. Several R/Bioconductor packages can be used to download and format the data [START_REF] Davis | GEOquery: a bridge between the Gene Expression Omnibus (GEO) and Bio-Conductor[END_REF][START_REF] Taminau | inSilicoDb: Access to the InSilico Database[END_REF]. We have chosen to encode in R our own functions. Our R script has been made available online, together with a user manual. Our formatting choices are described below.

In a GEO dataset two types of information are available for each sample. The first type consists of numeric values corresponding to a set of probes. The second type are character-type informations on the experimental setting. We have chosen to separate the two types into data matrices and information matrices.

In the data matrix, probes are associated to gene symbols with the use of different Bioconductor annotation packages according to the platform [START_REF] Carlson | Affymetrix Human Genome U133 Plus 2.0 Array annotation data[END_REF][START_REF] Carlson | Affymetrix Human Genome U133 Set annotation data (chip hgu133a)[END_REF][START_REF] Dunning | illuminaHumanv3.db: Illumina HumanHT12v3 annotation data (chip illumi-naHumanv3)[END_REF][START_REF] Dunning | illuminaHumanv4.db: Illumina HumanHT12v4 annotation data (chip illumi-naHumanv4)[END_REF]. After annotation, some symbols are duplicated. Several methods can be used to eliminate duplicates. We have chosen to keep the row with the largest interquartile range, as in [START_REF] Gentleman | genefilter: genefilter: methods for filtering genes from microarray experiments[END_REF], because we believe that this is the most statistically coherent choice. After annotation and reduction, the data matrix, with gene symbols as row names, and series numbers as column names, is saved as a single R object for future use. The information matrix has the same column names as the corresponding data matrix. Its rows correspond to the different fields.

Our merging function reduces data matrices to common row names. For information matrices, different sets of data usually have different information fields. This was taken into account when merging two information matrices, by indexing the rows of the merged matrix by the union of row names in the initial information matrices.

Two R/Bioconductor packages have recently been issued for merging GEO datasets [START_REF] Heider | virtualArray: a R/bioconductor package to merge raw data from different microarray platforms[END_REF][START_REF] Taminau | inSilicoMerging: Collection of Merging Techniques for Gene Expression Data[END_REF]. In [START_REF] Heider | virtualArray: a R/bioconductor package to merge raw data from different microarray platforms[END_REF], quantile discretization, normal discretization normalization, gene quantile normalization, median rank scores, quantile normalization (QN) are proposed. In [START_REF] Heider | virtualArray: a R/bioconductor package to merge raw data from different microarray platforms[END_REF][START_REF] Taminau | inSilicoMerging: Collection of Merging Techniques for Gene Expression Data[END_REF], the Batch Mean-Centering method, Distance-Weighted Discrimination, Z-score standardization, and the Cross-Platform Normalization method are proposed. An Empirical Bayes (EB) method is available in both packages. For the results reported here, only classical methods were used, and we consider them as sufficient to establish our main points, our focus being on overexpression or underexpression of genes, in different tissues or cancer types.

As in [START_REF] Tsodikov | Adjustments and measures of differential expression for microarray data[END_REF][START_REF] Warnat | Cross-platform analysis of cancer microarray data improves gene expression based classification of phenotypes[END_REF][START_REF] Breitling | Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments[END_REF], we have made the choice to use robust statistics [START_REF] Gibbons | Nonparametric statistical inference[END_REF][START_REF] Héritier | Robust methods in biostatistics[END_REF]. This implies changing the columns of a data matrix into distribution free values. The usually proposed transformation replaces the i-th value x i by its rank R i if x i is the R i -th smallest value in the column. However, ranks range between 1 and the number of rows. The problem is that different matrices may have different numbers of rows (gene symbols).

In order to get a unique range of values for all matrices, it seems preferable to use a scale free score. The simplest such score is the Empirical Cumulative Distribution Function (ECDF): its value at of p-values by the Benjamini-Yekutieli method [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] was systematically applied. Features were ranked from most to least significant, either by sorting p-values in increasing order, or by sorting the values of the test statistic instead. We considered as significant, any feature with a (FDR-corrected) p-value smaller than 5%.

x i is R i /n, if
Once a set of (significant) features had been selected, the corresponding rows were concatenated into single vectors. These vectors were taken as variables, and the samples as individuals, for a PCA. Figures 3 to 5

were obtained by projecting the samples as points onto the first principal plane, and differentiating their initial data matrices by colors. Precise R commands can be found in the user manual made available online.

Results

The 20 datasets that were downloaded from the GEO repository are detailed in Table 1. They were selected on a criterion of size (number of samples: 500 or more). The 20 matrices together amount to 17 745 samples.

To each study, a three-letter acronym was attached; these acronyms will be used in what follows.

acronym 1: Twenty GEO series have been chosen, coming from four different platforms. To each of them a three letters acronym was associated. The table gives the acronym, a recent reference, the GEO series number, the platform number. For the data matrix (or assayData), the number of symbols after annotation and reduction, and the number of columns (samples) are given. All 20 data matrices had 15 562 gene symbols in common.

In the results reported here, each data matrix has been transformed by replacing its column values, by the corresponding van der Waerden normal scores [7, p. 309]. Similar results were obtained when replacing column values by their empirical distribution function (see methods section).

The first treatment that was applied consisted in computing, for each dataset, the median of all rows, reduced to the 15 562 common gene symbols. This gave 20 vectors of length 15 562, the correlation matrix of which is given in Table 2. A positive (negative) correlation between vectors of size 15 562 is significant at threshold 5% if it is larger than 0.013 (smaller than -0.013); thus all correlations of Table 2 can be regarded as significant.

Figure 1 shows a factor analysis of the 20 variables. Fifteen of them can be clustered into four groups. • PMM, EPO, XMD, HBI, CCL, CCC. Among these six datasets, four are generalist studies involving different tissues and conditions (EPO, HBI, XMD, CCL); CCC concerns colon cancer, and PMM multiple myelomas. Observe that CCL, which was obtained under a platform different from the five others, has excellent correlations with them (between 0.74 and 0.92).

• WBS, PLE, HPS, HAV. All four correspond to blood RNA samples from healthy patients.

• MIL, AML, MDS. All three correspond to leukemias.

• GSC, MBI. These two matrices correspond to similar tissues (blood samples), and similar conditions (critical injuries and burn injuries). Moreover, they were produced on the same platform, by the same organization. Their excellent correlation (0.94) is not a surprise.

Three datasets, BEC, DLB, MMD have relatively good correlations with those of the above four groups (around 0.5), but no particular links with those groups, nor between themselves. The relative surprise comes from the weak correlations of PRS, and the negative correlations of PVA. Both come from blood RNA samples, and they could have been expected to be close to the WBS, PLE, HPS, HAV group. That PRS and PVA are far from any other matrix, can be explained by their inner heterogeneity. It is illustrated for PVA on Figure 2, where the values over features ALPP and CA4 are represented: samples separate into 4 clusters, according to over-or underexpression of the two genes. As an example, if PVA is split into samples for which the value of ALPP is positive (overexpression), or negative (underexpression), and the row medians are calculated over the two submatrices as before, a correlation of -0.69 is found: thus one half of PVA has a strong negative correlation with the other half. Similar results are obtained for many other features. We considered that the heterogeneity of PVA and PRS did not qualify them for merging.

For each matrix, we also computed all possible pairwise row correlations: 20 vectors of more than 121 millions of pair-correlations were obtained: this is the technique used to evaluate genes for cross-platform consistency of expression patterns in [START_REF] Parmigiani | A Cross-Study Comparison of Gene Expression Studies for the Molecular Classification of Lung Cancer[END_REF]. As expected, the correlation matrix had smaller values than that of Table 2. For instance, the correlation of CCL with XMD was 0.53 instead of 0.92, but still highly significant because of the large number of values.

Correlations between column medians or pair-correlations, is too crude a criterion to judge the homo- geneity of two datasets. As an example, GSC and MBI have an excellent median correlation of 0.94, and several good reasons to be similar. Yet, when each feature is tested for significant differences by the Kruskal-Wallis test, 14 800 significant features out of 18 260 are detected (see methods section for details). The same occurred for any pair of datasets: the distributions of rows had significantly different location parameters, for a majority of features. This means that, for a majority of genes, the ranks of their expressions in the first dataset are significantly smaller or larger than in the second. Since discrepancies appear to be observed between any two datasets, it must be decided whether they are due to actual biological information, or to a statistical artifact, induced by the experimental setting or the platform. For this, we focused on the dataset MIL (GSE13159 [START_REF] Haferlach | Clinical utility of microarray-based gene expression profiling in the diagnosis and subclassification of leukemia: report from the International Microarray Innovations in Leukemia Study Group[END_REF]), that has 2 096 samples. The samples were sorted into six submatrices, according to six keywords: Healthy (74 samples), ALL (acute lymphoblastic leukemia, 750 samples), AML (acute myeloid leukemia, 542 samples), CLL (chronic lymphocytic leukemia, 448 samples), CML (chronic myelogenous leukemia, 750 samples), MDS (myelodysplastic syndrome, 202 samples). Then the same treatments as before were applied. Firstly the six median columns were computed, and their correlation matrix was obtained (Table 3). The values are between 0.85 and 0.99, which is in the range of the best correlations of Table 2. As a control, we made a partition of the same matrix into 6 random subsets, with the same numbers of samples as above, and computed the correlation matrix in the same way. On the control random partition, all correlations were above 0.997. This proves that the partition into keywords does contain meaningful differences. Indeed, the samples projected as points onto the first principal plane, differentiating submatrices by colors (Figure 3). The two submatrices ALL (blue points) and CLL (brown points) are clearly separated from the rest.

-3 -2 -1 0 1 2 3 -3 -2 -
Healthy
Differences inside a given dataset can be induced by several factors. Two factors may not induce differences of the same order of magnitude. However, there is no statistical reason why a dataset like MIL should not be used as a whole, and many ways to verify that the observed differences correspond to actual biological information. Here is an example. Stirewalt et al. [START_REF] Stirewalt | identification of genes with abnormal expression changes in acute myeloid leukemia[END_REF] list a group of 7 genes displaying increased expression in acute myeloid leukemia samples: BIK, CCNA1, FUT4, IL3RA, HOMER3, JAG1, WT1. When a one-sided Wilcoxon test is applied to the submatrix AML versus the rest of MIL, those 7 genes are among the most significant: their p-values range between 6.7 × 10 -130 and 4.3 × 10 -42 . The most significant, HOMER3, ranks 54-th among the 20 184 features of MIL.

If observed differences between two datasets (like GSC and MBI) are of the same order of magnitude as differences inside a given dataset, such as caused by a significant factor (see figure 3), it can be admitted as statistically legitimate to merge the two datasets. That meaningful information can be obtained from the merging, remains to be proved. In the following experiments, matrices to be merged were selected in the clusters detected by factor analysis (Figure 1).

Our first experiment consisted in extracting samples corresponding to breast and lung tumors, from the three matrices CCL, EPO, and XMD. CCL has 56 samples of breast tumors, and 166 of lung tumors, EPO has 367 and 143, XMD has 32 and 152. Two matrices "Breast" and "Lung" were made by merging the six submatrices three by three, according to tissues. They had 18 466 features in common, by 455 samples for Breast, and 461 for Lung.

The Kruskal-Wallis test was run on the six separated submatrices, then on the two matrices Breast and Lung. The ten most significant symbols were extracted, and a PCA was run as before. The results are displayed on Figure 4. Significant symbols when the 6 matrices are separated (left panel) are different from significant symbols separating Breast and Lung (right panel). On the left panel, it is clear that the information on the dataset (CCL, EPO, or XMD) dominates the separation Breast vs. Lung: samples coming from CCL are on the left, from EPO on the right, from CCL in the middle. But on the right panel, - the two types of tumors are also clearly separated. Separators include GATA3 on the right side (Breast), IGF2BP3 on the left side (Lung). Two articles, among others, show the importance of GATA3 for breast cancer [START_REF] Zakaria | Identification of Estrogen-Related Genes in Breast Cancer: The Malaysian Context[END_REF][START_REF] Ma | The Cancer Genome Atlas: Clinical Applications for Breast Cancer[END_REF]. In [START_REF] Beljan | IMP3 can predict aggressive behaviour of lung adenocarcinoma[END_REF], the link of IGF2BP3 to lung cancer is explicitly stated. Further information was obtained by running a one-sided Wilcoxon test to detect symbols separating both types of tumor. Then the Molecular Signature database C2 [START_REF] Subramanian | Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles[END_REF] was searched for symbols matching them. Among the 20 genes found most significantly overexpressed in breast tumors by our test, 11 were inside genesets of C2 relative to breast cancers, and outside all genesets relative to lung tumors: EFHD1, IRX5, MUCL1, PRLR, PTGER3, RGL2, TRIL, TRPS1, VAV3, WWP1, ZG16B. Seven of these genes can be found in the G2SBC database [START_REF] Mosca | A multilevel data integration resource for breast cancer study[END_REF] and for 10 out of 11, we have found at least one reference relating it to breast cancer. Conversely, among the most significant genes for lung tumor, the following were found in C2 genesets related to lung and not in those related to breast: ALDH3B1, DARS, PRPSAP2, FAM96B, MBIP, LRRC20. The overexpression of ALDH3B1 in lung tumors has been reported in [START_REF] Marchitti | Aldehyde Dehydrogenase 3B1 (ALDH3B1): Immunohistochemical Tissue Distribution and Cellular-specific Localization in Normal and Cancerous Human Tissues[END_REF]. Santarius et al. [START_REF] Santarius | A census of overexpressed human cancer genes[END_REF] gives lists of genes, the overexpression of which is associated to different types of human cancers.

The genes detected as significantly overexpressed in Breast by our test, that were also among class III genes related to breast cancer in Table 1 of [START_REF] Santarius | A census of overexpressed human cancer genes[END_REF], were FGFR1, BAG4, MDM2, YWHAB, ZNF217. For Lung, they were EGFR, MET, YWHAZ, MYC, NKX2-1, DCUN1D1. These findings would require further confirmation over larger datasets. Yet they provide evidence that meaningful biological information can be extracted by merging generalist matrices such as CCL, EPO and XMD.

Our next experiment consisted in merging the two groups of blood RNA datasets, found to be homogeneous on the correlation analysis (Figure 1): HAV, HPS, PLE, WBS for healthy individuals, AML, MDS, for leukemias. The samples of MIL were separated into MILh (Healthy), and MILl (leukemias). The left panel of Figure 3 shows the first plane of the PCA for the same 22 features as in Figure 2, the 8 matrices being represented by different colors. It turns out that the samples corresponding to MILh are mixed on the representation with the other MIL points. Thus they were removed from the matrix "Healthy", whereas "Leukemia" was made by merging AML, MDS, and MILl. The Kruskal-Wallis test between Healthy and Leukemia, detected 16 977 significant features out of 17 691, among which 7 970 had a null p-value. The right panel of Figure 5 shows the PCA over 10 of them.

The one-sided Wilcoxon test was run to detect which symbols were significantly overexpressed in leukemias. For that test, a set of 4 191 symbols had a null p-value. A second set of symbol was extracted from C2: those appearing in leukemia-related genesets. The C2 set has 5 688 symbols, and the intersection with the first contains 1 617, which is highly significant for Fisher's hypergeometric test (P = 1.36 × 10 -51 ). The ten symbols found most significant for leukemia by the Wilcoxon test were RPL34, GABARAP, RPL36A, H2AFV, CSDE1, DNTTIP2, OPHN1, PABPC3, PNRC1, RPSA. Among those 10, 8 appeared in the leukemia-related genesets of C2. The symbol H2AFV is found in six of them. Another noteworthy result concerns the pair of genes NUP98-TOP1, shown to be related with leukemia in [START_REF] Gurevich | NUP98-Topoisomerase I acute myeloid leukemia-associated fusion gene has potent leukemogenic activities independent of an engineered catalytic site mutation[END_REF]. When genes are 10 ranked by decreasing order of significance, NUP98 and TOP1 have ranks 187 and 65 respectively, which confirms their link with leukemia.

Another experiment was run on the same matrices, by separating acute myeloid leukemia samples, from all other samples. Thus the same calculations that had been run inside MIL before, were repeated over a larger number of samples. The acute myeloid leukemia samples were taken from AML and MIL (1 076 samples), others were obtained by merging HAV, HPS, PLE, WBS, MDS, with the non-AML samples of MIL (6 010 samples). The one-sided Wilcoxon test of comparison was run. For the 7 genes signaled as overexpressed in AML by [START_REF] Stirewalt | identification of genes with abnormal expression changes in acute myeloid leukemia[END_REF], the results were much more significant as before: the least significant pvalue was that of BIK: 8.4 × 10 -35 , whereas FUT4 and HOMER3 had p-values below machine precision. Contrarily to the study that had been conducted inside MIL, a clear confirmation was also obtained for the genes reported by [START_REF] Stirewalt | identification of genes with abnormal expression changes in acute myeloid leukemia[END_REF] to be underexpressed in case of AML. Five of them were in the common features of our matrices, four had p-values smaller than 10 -100 for underexpression in AML. In particular, PELO and PLXNC1 who had not been found significantly underexpressed in the first experiment, now had p-values 3.5 × 10 -238 and 4.4 × 10 -168 in the test on merged matrices.

Discussion

A new set of R functions has been developed. Like other packages [START_REF] Davis | GEOquery: a bridge between the Gene Expression Omnibus (GEO) and Bio-Conductor[END_REF][START_REF] Taminau | inSilicoDb: Access to the InSilico Database[END_REF], it performs the usual formatting operations. It also offers new functionalities for sorting lists of datasets according to information keywords.

Various robust statistics techniques are encoded. The script and a user manual have been made available online. Using these R functions, a large scale study of 20 GEO datasets, totalling 17 745 samples, has been conducted.

Our first conclusion is that Haibe-Kains et al. [START_REF] Haibe-Kains | Inconsistency in large pharmacogenomic studies[END_REF] were right in observing that inconsistencies between datasets make it dangerous to merge them without precautions. The risk is to declare as biologically significant, observations which are actually statistical artifacts. The first precaution is to transform the data into distribution-free values, i.e. to use robust statistics. This implies replacing the data of each sample by their empirical distribution function, or some other distribution-free score [START_REF] Gibbons | Nonparametric statistical inference[END_REF][START_REF] Héritier | Robust methods in biostatistics[END_REF]. Even after data have been homogenized, important discrepancies remain. For this reason, checking comparability between studies before merging them is imperative. One possible measure of similarity (among others, see for instance [START_REF] Parmigiani | A Cross-Study Comparison of Gene Expression Studies for the Molecular Classification of Lung Cancer[END_REF]) for two datasets is the correlation between medians, which has been used here. Two sets of samples corresponding to different conditions inside one given homogeneous dataset usually have correlations of medians above 0.8 (see Table 3). Arguably, it can be considered that two different datasets can safely be merged, if all pairedcorrelations between medians are above 0.8. This is not always the case, even between datasets coming from the same tissues, obtained under the same platform (see Table 2). Further ways of investigating possible discrepancies involve multivariate statistics. Graphical methods include Factor Analysis, Principal Component Analysis, Discriminant Analysis [START_REF] Härdle | Applied Multivariate Statistical Analysis[END_REF]. Inference can be done using the robust equivalents of usual normal-sample methods, i.e. Wilcoxon test instead of Student's t-test, Kruskal-Wallis instead of one-way anova, etc. When repeatedly applying such a test to a set of symbols, a False Detection Rate (FDR) correction must be applied to the p-values. We have chosen the Benjamini-Yekutieli method [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. Our observation was that, even after FDR correction, the tests usually detect a sizeable proportion of all symbols as significant for discrimination, either between several different datasets, or between different types of samples within the same dataset. We believe that relevant biological information can be obtained from applying a discriminating test, then ranking features according to their degree of significance, i.e. ordering the values obtained over each feature by the test statistic. In the cases considered here (breast tumors against lung tumors, healthy blood samples against leukemias, acute myeloid leukemia against other blood RNA samples), it was observed that among the most significant symbols, a large proportion of them were already known as being related to the corresponding cancers. This can be viewed as evidence that meaningful biological information can be extracted by merging different datasets. We believe that important new findings could be obtained by the same method, being aware that a statistical listing of significant symbols does not necessarily imply that all listed symbols correspond to true biological information. Such a list must necessarily be expert-curated for biochemical validation.

Figure 1 :

 1 Figure 1: Factor analysis of median columns for 20 datasets. The 20 variables are projected onto the first principal plane of the PCA. Four clusters are identified.

Figure 2 :

 2 Figure 2: Values of PVA on CA4 versus ALPP. Samples separate into 4 clusters, according to over-or underexpression of the two genes.

Figure 3 :

 3 Figure 3: Dataset MIL, partitioned into 6 submatrices according to keywords Healthy, ALL, AML, CLL, CML, MDS. PCA of the 22 symbols with Kruskal-Wallis p-value under 10 -300 : SOX4, SYNGR2, ER-LIN1, FAH, C7orf23, PSMA6, RTN3, UHRF1, ADAM28, BLK, FUCA2, CD79A, ADA, MYL6B, HEBP1, LEF1-AS1, LEF1, AFF3, COL9A2, MICALL2, MPO, PPM1K. Samples are represented by points, with six different colors.

Figure 4 :

 4 Figure 4: Principal component analysis of two assays of 10 symbols in 6 submatrices, extracted from CCL, EPO, and XMD according to keywords "Breast" and "Lung". The six submatrices are separated on the right panel, they have been merged on the right panel. In each case the 10 most significant features for the Kruskal-Wallis test are taken as variables. The two sets of 10 symbols are disjoint. Samples are represented by red points (Breast) or blue points (Lung).

Figure 5 :

 5 Figure 5: On the left panel, PCA of HAV, HPS, PLE, WBS, MILh (Healthy = red points), and AML, MDS, MILl (Leukemia = blue points), for the same 22 features as in Figure 2. On the right panel, HAV, HPS, ML, WBS, have been grouped into Healthy, AML, MDS, MILl into Leukemia. The assay is made of 10 random features among the 7 970 having a null Kruskal-Wallis p-value.

Table

		reference series number platform number symbols (rows) samples (columns)
	EPO	[22]	2109	570	20 184	2158
	PMM	[23]	2658	570	20 184	559
	AML	[24]	6891	570	20 184	537
	HBI	[25]	7307	570	20 184	677
	MIL	[26]	13159	570	20 184	2 096
	MDS	[27]	15061	570	20 184	870
	PLE	[28]	20142	6947	19 626	1 240
	MMD	[29]	24080	570	20 184	559
	DLB	[30]	31312	570	20 184	498
	PRS	[31]	33828	10558	20 768	881
	CCL	[32]	36133	15308	18 722	917
	BEC	[33]	36192	6947	19 628	911
	WBS	[34]	36382	6947	19 628	991
	GSC	[35]	36809	570	18 260	812
	MBI	[36]	37069	570	18 260	590
	CCC	[37]	39582	570	20 184	566
	PVA	[38]	48152	6947	19 628	705
	HPS	[39]	48348	6947	19 628	734
	XMD	[40]	48433	570	20 184	823
	HAV	[41]	48762	6947	19 628	621

Table 2 :

 2 For each of the 20 data matrices of Table1, the median column value of each gene symbol was computed. This gave 20 vectors with length 15 562 (number of common symbols). The table gives pairwise correlations between the 20 vectors.

Table 3 :

 3 The data matrix MIL was partitioned according to the 6 keywords Healthy, ALL, AML, CLL, CML, MDS. For each of the six submatrices, the median column of each feature was computed. This gave 6 vectors with length 20 184 (number of symbols in MIL). The table gives the correlations of the 6 vectors. these differences were detected by the Kruskal-Wallis test: out of the 20 184 features, 18 301 were found significant. Twenty-two features had Kruskal-Wallis p-value below 10 -300 : SOX4, SYNGR2, ERLIN1, FAH,

			ALL AML CLL CML MDS
	Healthy	1.00	0.92	0.96	0.86	0.98	0.99
	ALL	0.92	1.00	0.95	0.91	0.90	0.91
	AML	0.96	0.95	1.00	0.89	0.96	0.97
	CLL	0.86	0.91	0.89	1.00	0.84	0.85
	CML	0.98	0.90	0.96	0.84	1.00	0.98
	MDS	0.99	0.91	0.97	0.85	0.98	1.00
	C7orf23, PSMA6, RTN3, UHRF1, ADAM28, BLK, FUCA2, CD79A, ADA, MYL6B, HEBP1, LEF1-AS1,
	LEF1, AFF3, COL9A2, MICALL2, MPO, PPM1K. A PCA of the corresponding rows of MIL was run, and

a pdf file sagd manual.pdf: user manual for the R functions.

Acknowledgements

BY, KC, and JJF acknowledge financial support from Laboratoire d'Excellence TOUCAN (Toulouse Cancer).

Additional Files

Additional material has been provided as a compressed directory available online: http://ljk.imag.fr/membres/Bernard.Ycart/publis/sagd.tgz It contains:

1. a R script file sagd.r: the R functions implementing the method described here,