BACKGROUND
a crucial issue in genomic studies is the elucidation of how genes change expression and interact as a consequence of external/internal stimuli such as an illness, drug administration, hormone stimuli, etc.
microarray technology makes it possible to monitor simultaneously a large number of gene transcripts through a series of different experimental conditions.
in particular, microarray time series studies are essential to understand the dynamics of biological events at the molecular level.
a first necessary step in order to limit the analysis to those genes that change expression over time is to select differentially expressed transcripts.
selection methods proposed in the literature usually deal with the comparison of static  rather than dynamic conditions, and are based on statistical tests  <cit> .
these methods test the significance of the differential expression gene by gene.
at least two replicates for each of the conditions to be tested are necessary, but a higher number is required to have reliable results.
in time series experiments, in which gene expression is monitored over time, it is necessary to test differential expression at different sampling times.
anova or anova based procedures  <cit>  have been proposed to this purpose.
however, since in time series experiments replicates are often available only for a limited number of samples, anova tests are seldom applicable.
for this reason, differentially expressed genes in time series experiments are often selected using an empirical constant fold change threshold  <cit> .
this is far from ideal, since it is based on an arbitrary choice , which does not take into account the characteristics of the measurement error.
when the number of the replicates is not sufficient to apply traditional statistical tests, alternative methods need to be applied.
two methods based on a fit of the time series were recently proposed in the literature  <cit> .
these methods fit the time series expression profiles using respectively polynomials and splines.
comparison between time series is based respectively on model parameters and goodness of fit.
both methods are really general and do not require any replicates; however, it is not clear the role of the number of available samples on their performance.
here we propose methods  <dig> and  <dig> able to select differentially expressed gene profiles in data-poor conditions, based on a model of the experimental error.
their performance is investigated in comparison to method  <cit>  , based on splines fit, using synthetic time series of different length.
finally, a case study on insulin treated muscle cells is presented to better appreciate the implementation aspects of methods  <dig> and  <dig>
methods
selection strategy
let's call xt and xc the log-expression measurements in treated  and control  cultures, available for a generic gene x at time sample tk .
log expression measurement are used, as in  <cit> , because the signal is considered proportional to the log of the measurements, the error is considered log-additive, and the large range of expression intensities makes the log-expression practical.
the rationale adopted to label a gene x as differentially expressed in condition t vs c is described in details for methods  <dig> and  <dig> and is briefly reviewed for method  <dig>  since we refer to  <cit>  for further details.
method 1
the deviation of expression of gene x in t and c is calculated for each sample tk as:
d = xt - xc
the gene is considered differentially expressed in t vs c if |d| exceeds a threshold θd in at least one sampling time tk :
|d| >θd
where θd is determined in correspondence to a significance level α, based on the null hypothesis distribution of d. this distribution is modeled from d values calculated by equation  <dig>  with xt and xc measured on experimental replicates , which provide a situation where genes are not differentially expressed, so that the null hypothesis is verified.
therefore, replicates of at least one time sample are necessary, to apply method  <dig>
method 2
the area a bounded by the two expression profiles t and c  is calculated for each gene x as the sum of the contributions of partial areas from consecutive pairs of samples:
a=∑k=1m−1ak
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqgbbqqcqgh9aqpdaaewbqaaiabbgeabnaabaaaleaacqqgrbwaaeqaaaqaaiabbugarjabg2da9iabigdaxaqaaiabb2eanjabgkhitiabigdaxaqdcqghris5aogaaczcaiaaxmaadaqadaqaaiabiodazagaayjkaiaawmcaaaaa@3da0@
each contribution ak is calculated from the deviation of expression in t and c , as:
if  sign)=sign)then  ak=+d)⋅2else  ak=⋅dd+d+d⋅dd+d)⋅2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqaaewabaaabagaeeysakkaeeoraykaeeiiaaiaeeiiaaiaee4camnaeeyaakmaee4zacmaeeoba4maeiikagiaeeizaqmaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqggpaqkcqggpaqkcqgh9aqpcqqgzbwccqqgpbqacqqgnbwzcqqgubgbcqggoaakcqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapiabcmcapaqaaiabbsfaujabbieaijabbweafjabb6eaojabbccagiabbccagiabbgeabnaabaaaleaacqqgrbwaaeqaaogaeyypa0zaasaaaeaacqggoaakcqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasmaey4kasiaegymaedabeaakiabcmcapiabgucariabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaeiykakiaeyyxictaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqghsislcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabikdayaaaaeaacqqgfbqrcqqgmbatcqqgtbwucqqgfbqrcqqggaaicqqggaaicqqgbbqqdawgaawcbagaee4aasgabeaakiabg2da9maalaaabawaaewaaeaacqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasmaey4kasiaegymaedabeaakiabcmcapiabgwsixpaalaaabagaeeizaqmaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqggpaqkaeaacqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasmaey4kasiaegymaedabeaakiabcmcapiabgucariabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakcaaiabgucariabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaeyyxic9aasaaaeaacqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwacqghrawkcqaixaqmaeqaaogaeiykakiaey4kasiaeeizaqmaeiikagiaeeidaq3aasbaasqaaiabbugarbqabagccqggpaqkaaaacagloagaayzkaagaeyyxictaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqghsislcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabikdayaaaaagaaczcaiaaxmaadaqadaqaaiabisda0agaayjkaiaawmcaaaaa@d0ee@
the gene x is considered differentially expressed in t vs c if the following inequality holds:
a >θa
where θa is a threshold to be determined, in correspondence to a significance level α, based on the null hypothesis distribution of a, i.e.
the distribution of a derived from experimental replicates .
method 3
for each gene x, the expression profiles c and t are fitted using natural cubic splines.
the null hypothesis is that both c and t time series share the same fit, the alternative hypothesis is that the fits are not equal.
as a first step of the method, the same spline function is fitted simultaneously to the two profiles c and t and the sum of squares residuals ss <dig> is calculated; then two different spline functions are fitted separately to time series c and t and the sum of squares residuals ss <dig> is computed.
to assess differentially expressed genes, the goodness of model fit under the null hypothesis is compared to that under the alternative hypothesis, by calculating a statistic f as:
f=ss0−ss1ss0
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqggbgrcqgh9aqpdawcaaqaaiabbofatjabbofatnaacaaaleqabagaegimaadaaogaeyoei0iaee4uamlaee4uam1aawbaasqabeaacqaixaqmaaaakeaacqqgtbwucqqgtbwudaahaawcbeqaaiabicdawaaaaagccawljagaaczcamaabmaabagaegonaydacagloagaayzkaaaaaa@3e08@
gene x is considered differentially expressed if the following inequality holds:
f >θf
where θf is a threshold to be determined, in correspondence to a significance level α, based on the null hypothesis distribution of f.
null hypothesis distribution model
both methods  <dig> and  <dig> are based on a threshold derived from the null hypothesis distribution of d and a, respectively, obtained using available replicates.
we denote as "available replicates" replicates available for a subset of time samples .
let's suppose two experimental replicates are available for a generic time tk and a generic experimental condition .
by assuming a log-additive error model as in  <cit> , the log-expression measurement of each gene x in replicates a and b, can be expressed as:
{xa=μ+εaxb=μ+εb
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagabeqaauaabeqaceaaaeaacqqg4baedawgaawcbagaeeyyaegabeaakiabg2da9ggaciab=x7atjabgucariab=v7alnaabaaaleaacqqghbqyaeqaaagcbagaeeieag3aasbaasqaaiabbkgaibqabagccqgh9aqpcqwf8oqbcqghrawkcqwf1oqzdawgaawcbagaeeoyaigabeaaaaaakiaawuhaaiaaxmaacawljawaaewaaeaacqai4aaoaiaawicacaglpaaaaaa@450f@
where μ represents the actual gene expression  and εa, εb are two realizations of the error variable ε, monitoring both technical and biological variability.
the indices for condition and time tk are omitted here because we refer to pair of replicates available for a generic time sample tk.
null hypothesis distribution of variable d
to quantify the deviation of expression d under the null hypothesis, equation  <dig> is applied to available replicates as:
dh <dig> = xa - xb = εa - εb
different distribution models  are used to fit the set of dh <dig> values obtained by applying equation  <dig> to all genes and available replicates.
the best model is selected based on the goodness of fit and the parameters precision, and is used as the null hypothesis distribution of d to determine θd to be used in equation  <dig>
to determine the threshold in correspondence to a significance level α, method  <dig> uses a model to fit the observed statistics rather than using quantiles.
the reason for this choice is that the lack of a sufficient number of observations from available replicates renders the determination of appropriate thresholds difficult when low significance levels are chosen, as often the case in microarray studies.
if a sufficient number of replicates is available to guarantee a good threshold setting at the desired significance level α, it may be preferable to use quantiles.
null hypothesis distribution of variable a
at least two replicates for each time sample would be necessary to derive a distribution under the null hypothesis from the data.
since we address selection of differentially expressed genes in data-poor condition, i.e.
a sufficient number of replicates is not available, a monte carlo procedure is used to derive the null distribution of a. first, dh <dig> distribution is derived from d values obtained from available replicates as described above.
then, b profiles of length m are sampled from dh <dig>  under the hypothesis that the error at different time samples is independent and identically distributed.
subsequently, b values of ah <dig> are calculated from these profiles.
finally, different distribution models  are used to fit the entire set of ah <dig> values and the best model is chosen based on goodness of fit and parameter precision.
this model is used as the null hypothesis distribution of a to determine θa to be used in equation  <dig>
as for method  <dig>  a model is fitted to the observed statistics rather than using quantiles to determine the threshold in correspondence to a significance level α.
null hypothesis distribution of variable f
the null distribution of f is obtained using bootstrap.
see  <cit>  for details.
intensity dependency of error
in affymetrix chips it is well known that dh <dig> has an intensity dependent distribution  <cit> .
in particular, analysis of technical replicates of affymetrix human chip has shown that the standardized variable sh <dig> :
sh0=dh0sddh0=xa−xbsddh0;x¯=xa+xb2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqabeqacaaabagaee4cam3aawbaasqabeaacqqgibascqaiwaamaagccqgh9aqpdawcaaqaaiabbsgaknaacaaaleqabagaeeisagkaegimaadaaagcbagaee4uamlaeeiraq0aasbaasqaaiabbsgakjabbieaijabbcdawaqabagccqggoaakcuqg4baegaqeaiabcmcapaaacqgh9aqpdawcaaqaaiabbiha4naabaaaleaacqqghbqyaeqaaogaeyoei0iaeeieag3aasbaasqaaiabbkgaibqabaaakeaacqqgtbwucqqgebardawgaawcbagaeeizaqmaeeisagkaegimaadabeaakiabcicaoiqbbiha4zaaragaeiykakcaaiabcuda7aqaaiqbbiha4zaaragaeyypa0zaasaaaeaacqqg4baedawgaawcbagaeeyyaegabeaakiabgucariabbiha4naabaaaleaacqqgibgyaeqaaagcbagaegomaidaaaaacawljagaaczcamaabmaabagaegymaejaegimaadacagloagaayzkaaaaaa@5f0b@
has an intensity independent distribution  <cit> .
therefore, in case of data showing intensity dependency of the variable dh <dig>  it is convenient to model sh <dig> distribution, as indicated in the additional file  <dig>  to derive the threshold θd to be used in equation  <dig>  consistently, the values of d observed from the data are standardized before applying equation  <dig>  if method  <dig> is used:
s=dsddh0);x¯=xt+xc2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqabeqacaaabagaee4camnaeiikagiaeeidaq3aasbaasqaaiabbugarbqabagccqggpaqkcqgh9aqpdawcaaqaaiabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakcabagaee4uamlaeeiraq0aasbaasqaaiabbsgakjabbieaijabicdawaqabagccqggoaakcuqg4baegaqeaiabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaeiykakcaaiabcuda7aqaaiqbbiha4zaaragaeiikagiaeeidaq3aasbaasqaaiabbugarbqabagccqggpaqkaagaeyypa0zaasaaaeaacqqg4baedaahaawcbeqaaiabbsfaubaakiabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaey4kasiaeeieag3aawbaasqabeaacqqgdbwqaagccqggoaakcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabikdayaaacawljagaaczcamaabmaabagaegymaejaegymaedacagloagaayzkaaaaaa@6586@
analogously, if method  <dig> is used on data showing intensity dependency of the variable dh <dig>  ah <dig> and θa are derived using sh <dig> and the values of a observed from the data  are calculated using s instead of d.
threshold setting
once the null hypothesis distribution of d, a and f are obtained, thresholds θd, θa and θf are determined in correspondence to a significance level α. rather than fixing it a priori, α can be optimized based on a variety of criteria aiming to control the family wise error rate   <cit> , or the false discovery rate   <cit>  or a compromise between false positive and false negative classification  <cit> .
as an example, let's focus on a criterion based on the control of fdr, defined as the expected proportion of false positive classification  among the number sα of genes selected as differentially expressed, using significance level α:
fdr=e
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqggbgrcqqgebarcqqgsbgucqgh9aqpcqqgfbqrdawadaqaamaalaaabagaeeoraykaeeiuaafabagaee4uam1aasbaasqaaggaciab=f7ahbqabaaaaagccaglbbgaayzxaagaaczcaiaaxmaadaqadaqaaiabigdaxiabikdayagaayjkaiaawmcaaaaa@3e04@
in case of numerous sets as for microarrays, fdr is well approximated by
fdr≈esα
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqggbgrcqqgebarcqqgsbgucqghijyudawcaaqaaiabbweafjabcufabjabdaeagjabdcfaqjabc2fadbqaaiabdofatnaabaaaleaaiigacqwfxoqyaeqaaaaakiaaxmaacawljawaaewaaeaacqaixaqmcqaizawmaiaawicacaglpaaaaaa@3f45@
calculating fdr requires the estimate of the expected number of false positives, obtained as the product of α by the number n <dig> of non differentially expressed genes:
e  = n <dig> · α
n <dig> is unknown and is estimated using the bootstrap procedure described in  <cit> .
fdr is calculated for a range of significance levels α and the significance level that guarantees the desired fdr is then used to select differentially expressed genes.
since method  <dig> applies m tests  for each gene, the significance level α for method  <dig> is corrected by applying Šidák correction  <cit>  in order to account for multiple testing.
simulation
three different experimental conditions with a number of time samples m =  <dig>   <dig>   <dig> were simulated.
<dig> synthetic data sets were generated for each experimental condition, each consisting of  <dig> profiles:  <dig> simulated as differentially expressed and the remaining as random noise.
in both cases, the deviation of expression in t vs c was generated at each sampling time tk as standardized deviation:
s ~ n ∀k =  <dig>  ..., m
where σ <dig> was set equal to  <dig> and μk =  <dig>  for not differentially expressed genes; while, for differentially expressed genes, plausible profiles were obtained by modeling μk as dependent on μk- <dig> according to a first order markov model , with the only constraint of being greater than  <dig>  for at least one time samples.
samples k =  <dig> were generated twice for each gene, so as to provide replicates useful to apply method  <dig> and  <dig>  these replicates were included also in the analysis for method  <dig>  simulated data were used to test the performance of the methods in different experimental conditions.
after the null hypothesis distributions of variables d, a and f are modeled as described above, a significance level α had to be fixed to determine the confidence thresholds θd, θa and θf to be used in equations  <dig>   <dig> and  <dig> respectively.
we compared the performance of the three methods across the entire range of significance level α by using precision  vs recall , and curves of observed false positives divided by the number of selected genes  vs number of selected genes.
moreover, we compared the average results obtained with the  <dig> methods by setting α in correspondence to a desired fdr of  <dig> , in terms of number of selected genes and observed false positives divided by the number of selected genes.
all measurements were averaged across the  <dig> simulations.
insulin case study
to better appreciate some characteristics of methods  <dig> and  <dig> related to the experimental error modeling, the analysis of the null hypothesis distribution of the variables d and a  was applied to a real case study on rat muscle cells treated with insulin.
the study was performed in vitro, on muscle l <dig> rat cell line.
cells were treated with insulin at time 0+, just after the collection of a first baseline sample at time 0; eight samples were harvested every hour during eight hours insulin stimulation.
a control experiment was also performed in order to be able to distinguish between insulin effect and biological processes of different nature, which take place in the culture simultaneously to insulin induced processes .
a total of twenty affymetrix chips rg_u34a  were hybridized using four replicates of the basal sample, eight samples collected from the control culture, and eight samples collected from the treated culture.
standard affymetrix mas  <dig> software  <cit>  was used for data pre-processing.
RESULTS
simulation
an example of simulated data with m =  <dig> is shown in figure  <dig>  as clusters obtained by using self organizing maps  <cit>  and pearson correlation as similarity measure on one of the simulated data-set .
the average profile is shown for each cluster, together with standard deviation bars; the number of genes for each cluster is also reported.
a variety of profiles are represented, such as genes differentially expressed in one or few peaks apparently uncorrelated , profiles that show characteristic bumps and waves of different length  or consistent trends along the time series, accompanied by few  or numerous  time samples with average absolute value greater than the error standard deviation σ .
results obtained by applying methods  <dig>   <dig> and  <dig> to  <dig> simulated data sets are shown as average precision at different ranges of recall intensities in figure  <dig>  left panels.
average and standard deviations across the  <dig> simulations are reported.
method  <dig> slightly outperforms method  <dig> for precision higher than  <dig> , since it provides higher values of both precision and recall in all the different sampling conditions.
for values of precision lower than  <dig> , the two curves are almost superimposable.
both methods  <dig> and  <dig> work better than method  <dig> for short time series , while for long time series  the performance of the methods  <dig> and  <dig> is similar, with method  <dig> slightly superior to methods  <dig>  the areas under the curves  are also reported in figure  <dig> .
to better appreciate the ability of the three methods to select differentially expressed genes at low fdr, figure  <dig>  shows the number of observed false positives divided by the number of selected genes  for different numbers of selected genes in the range 0–300: method  <dig> provides number of false positives divided by the number of selected genes always lower than that obtained with method 1; method  <dig> is slightly superior to method  <dig> for m =  <dig>  but inferior to both methods  <dig> and  <dig> for m =  <dig> and m =  <dig>  the performance of the three methods with respect to each other is not affected by the number of differentially expressed genes .
insulin case study
error
in this section we present the estimate of the null hypothesis distribution of variable d  and variable a .
figure  <dig> shows the values of dh <dig>   and its standardized value sh <dig>  plotted vs the average intensity of gene expression values .
figure  <dig>  shows the average variance of dh <dig> and sh <dig> vs the intensity range of gene expression discretized in intervals of constant size.
results confirm that while dh <dig> has an intensity-dependent distribution, sh <dig> has not  <cit>  .
the entire set of sh <dig> values was fitted using different distribution models .
the gaussian mixture model with n =  <dig> was chosen as the best model.
model parameters and their precision are shown in table  <dig>  monte carlo procedure was then used to derive the null distribution of a: b profiles of length m =  <dig> were sampled from sh <dig> distribution and b =  <dig> values of ah <dig> were calculated from these profiles.
three different distribution models  were used to fit the entire set of ah <dig> values.
the best fit for ah <dig> was obtained with gamma distribution .
gene selection
we applied the three methods to the data, adopting the false discovery rate as criterion for threshold setting.
table  <dig> shows the number of probe-sets selected by each of the three methods for fdr equal to  <dig> ,  <dig> ,  <dig> ,  <dig>  respectively.
results confirm those obtained by using simulation, i.e.
the ability of method  <dig> to select a higher number of differentially expressed genes at controlled fdr for short time series.
since focus here is methodology, biological results are not discussed further; confirmation studies and biological interpretation will be presented in a different article.
discussion
in this work we evaluated the performance of two selection methods here proposed to be applied in time series studies in data-poor conditions, i.e.
when the number of available replicates does not make possible or practical the use of standard statistical methods.
we also tested the two methods in comparison with a third method from the literature.
method  <dig> compares samples time by time using a statistically based fold change threshold derived from a null hypothesis distribution of variable d. to this purpose, replicates of at least one time sample are necessary.
since the threshold is derived based on the experimental variability, method  <dig> accounts for the error characteristics, e.g.
its intensity dependence; therefore, for example using affymetrix chips, genes expressed at high intensities are not penalized with respect to genes at low intensities, which show a higher variability .
although method  <dig> improves upon the use of a constant empirical fold change threshold, it considers time samples independently to each other, which is not a realistic assumption in time series studies.
method  <dig> calculates the area of the region bounded by the time series expression profiles to be compared and considers the gene differentially expressed if this area exceeds a threshold based on a model of the experimental error; therefore, besides accounting for the error distribution, it considers the entire expression profile and not single time samples.
also this method needs replicates for at least one time sample in order to derive the null hypothesis distribution.
both methods  <dig> and  <dig> assumes as working hypothesis that the error at different time points is independent and identically distributed.
this hypothesis, on our experience, is usually verified on real data .
however, for some experimental settings there may be a dependency of the error on time.
in this latter case, it would be more appropriate to perform replicates at different time samples, covering the duration of the experiment, and use method  <dig> with time dependent threshold settings.
methods  <dig> and  <dig> can be applied to compare time series from  <dig> different experimental conditions or a time series vs its baseline , by defining the deviation of expression of gene × for each sample tk  as the deviation between treated at time tk and the baseline.
moreover, if the sampling grid is different for the two time series, method  <dig> can be easily generalized by generating ah <dig> distribution from b time series with appropriate sampling grid.
method  <dig>  <cit>  uses natural cubic splines to fit time series expression data using a null hypothesis and an alternative hypothesis model, and performs a statistical test on the sum of squares residuals obtained using the two models to assess differential expression.
it therefore, besides considering the entire expression profiles, implicitly considers time dependencies.
moreover, it does not need any experimental replicate, which makes it practical for a wide range of time series microarray experiments.
we tested the performance of the three methods using synthetic data to assess their validity over a range of experimental conditions, specifically the length of the analyzed time series.
in the simulation, we used a markov model to generate plausible data, accounting for the dependencies of time samples and not biased toward one of the methods.
taking for example μk equal to an arbitrary constant for a random subset of time samples tk  in equation  <dig>  would not have accounted for time dependencies, and would have generated non realistic oscillation in the profiles, thus penalizing method  <dig> which is based on a model fit, with respect to the other methods, in particular method  <dig>  which applies a threshold on each time sample.
moreover, simulated profiles  represent a variety of possible situations in time series expression, such as profiles characterized by one or few peaks, waves of different length or consistent trends along the time series.
results on simulated data showed that method  <dig> outperforms method  <dig> independently from the length of the time series being analyzed, probably because, as method  <dig> is based on single sample comparisons, it is particularly sensitive to random fluctuations due to the noise, thus resulting in a larger number of false positives.
method  <dig> constitutes an improvement with respect to method  <dig> since the entire expression profile is considered simultaneously and this allows better distinguishing between consistent differences in expression profiles and random oscillations, thus resulting in a lower false negative rate.
method  <dig>  as method  <dig>  considers the entire expression profile, but performs better than method  <dig>  only for long time series .
looking at the ability of the methods to classify profiles with particular characteristics, we observed that method  <dig> works better than method  <dig> to detect differentially expressed genes that show just one or two peaks as differentially expressed.
on the opposite, method  <dig> works better than method  <dig> in detecting profiles characterized by waves of length greater than three samples, so as profiles that show a characteristic increasing/decreasing trend.
method  <dig>  as method  <dig>  is better in detecting bumps and consistent trends in the profiles, than in detecting isolated peaks.
however, it needs long time series  to perform better than method  <dig>  probably because the fit are more reliable when performed on long time series; method  <dig> was in fact proposed by the authors on real case studies with more than  <dig> samples.
these results are certainly of interest to address the choice of the algorithm to be used in data-poor time series expression studies, depending on the availability of replicates and on the length of the time series.
CONCLUSIONS
microarray time series studies are essential to understand the dynamics of biological molecular events.
in order to limit the analysis to those genes that change expression over time, it is necessary to select differentially expressed transcripts.
due to the high cost of microarrays, experiments are often performed without replication; therefore, traditional statistical methods can't be applied.
here we evaluate the performance of two selection methods applicable in data poor conditions, based on: a statistically based threshold on individual samples; a statistically based threshold to be applied on the area of the region bounded by the time series expression profiles to be compared.
application on a real data set on insulin regulation on muscle cells, obtained using affymetrix chips, revealed as the error analysis performed using methods  <dig> and  <dig> may be useful to detect error characteristics such as intensity dependencies and to properly address these feature by standardization.
we evaluated methods  <dig> and  <dig> performance using simulated data with a different number of available samples and compared these performance with those obtained using method  <dig>  based on a splines fit of time series profiles.
the results outlines that the two error based methods  <dig> and  <dig> work better than method  <dig> with short time series experiments, while method  <dig> works better than methods  <dig> and  <dig> with long time series experiments.
these results might help to optimize the choice of the algorithm to be used in different experimental conditions.
a preliminary version of method  <dig>  implemented in r, is available at  <cit> .
authors' contributions
bdc conceived the study and performed data analysis under the guidance and supervision of gt.
skn was responsible for the overall coordination of microarray experiments realization.
ljg performed microarray experiments.
cc was responsible for the overall project coordination.
all authors read and approved the final manuscript.
supplementary material
additional file 1
the .pdf document contains a detailed description of the procedure used to standardize the variable dh <dig> so as to account for the error intensity dependent distribution, and a detailed description of the procedure used to simulate the data, based on a first order markov model.
click here for file
acknowledgements
this study was supported by ministero dell'università e della ricerca scientificae tecnologica , and by national institute of health, grant roidk <dig>
this article has been published as part of bmc bioinformatics volume  <dig>  supplement  <dig>  2007: italian society of bioinformatics : annual meeting  <dig>  the full contents of the supplement are available online at .
figures and tables
all measurements were averaged across the  <dig> simulations; standard deviations are reported in parenthesis.
number of selected genes and number of false positives divided by the number of selected genes, obtained by applying the  <dig> methods on simulated data .
number of genes selected applying the  <dig> methods on real data .

 a crucial issue in genomic studies is the elucidation of how genes change expression and interact as a consequence of external/internal stimuli such as an illness, drug administration, hormone stimuli, etc. microarray technology makes it possible to monitor simultaneously a large number of gene transcripts through a series of different experimental conditions. in particular, microarray time series studies are essential to understand the dynamics of biological events at the molecular level. a first necessary step in order to limit the analysis to those genes that change expression over time is to select differentially expressed transcripts. selection methods proposed in the literature usually deal with the comparison of static rather than dynamic conditions, and are based on statistical tests. these methods test the significance of the differential expression gene by gene. at least two replicates for each of the conditions to be tested are necessary, but a higher number is required to have reliable results. in time series experiments, in which gene expression is monitored over time, it is necessary to test differential expression at different sampling times. anova or anova based procedures have been proposed to this purpose. however, since in time series experiments replicates are often available only for a limited number of samples, anova tests are seldom applicable. for this reason, differentially expressed genes in time series experiments are often selected using an empirical constant fold change threshold. this is far from ideal, since it is based on an arbitrary choice, which does not take into account the characteristics of the measurement error. when the number of the replicates is not sufficient to apply traditional statistical tests, alternative methods need to be applied. two methods based on a fit of the time series were recently proposed in the literature. these methods fit the time series expression profiles using respectively polynomials and splines. comparison between time series is based respectively on model parameters and goodness of fit. both methods are really general and do not require any replicates; however, it is not clear the role of the number of available samples on their performance. here we propose methods and able to select differentially expressed gene profiles in data-poor conditions, based on a model of the experimental error. their performance is investigated in comparison to method, based on splines fit, using synthetic time series of different length. finally, a case study on insulin treated muscle cells is presented to better appreciate the implementation aspects of methods and methods selection strategy let's call xt and xc the log-expression measurements in treated and control cultures, available for a generic gene x at time sample tk. log expression measurement are used, as in, because the signal is considered proportional to the log of the measurements, the error is considered log-additive, and the large range of expression intensities makes the log-expression practical. the rationale adopted to label a gene x as differentially expressed in condition t vs c is described in details for methods and and is briefly reviewed for method since we refer to for further details. method 1 the deviation of expression of gene x in t and c is calculated for each sample tk as: d = xt - xc     the gene is considered differentially expressed in t vs c if |d| exceeds a threshold θd in at least one sampling time tk: |d| >θd     where θd is determined in correspondence to a significance level α, based on the null hypothesis distribution of d. this distribution is modeled from d values calculated by equation with xt and xc measured on experimental replicates, which provide a situation where genes are not differentially expressed, so that the null hypothesis is verified. therefore, replicates of at least one time sample are necessary, to apply method method 2 the area a bounded by the two expression profiles t and c is calculated for each gene x as the sum of the contributions of partial areas from consecutive pairs of samples: a=∑k=1m−1ak      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqgbbqqcqgh9aqpdaaewbqaaiabbgeabnaabaaaleaacqqgrbwaaeqaaaqaaiabbugarjabg2da9iabigdaxaqaaiabb2eanjabgkhitiabigdaxaqdcqghris5aogaaczcaiaaxmaadaqadaqaaiabiodazagaayjkaiaawmcaaaaa@3da0@ each contribution ak is calculated from the deviation of expression in t and c, as: if  sign)=sign)then  ak=+d)⋅2else  ak=⋅dd+d+d⋅dd+d)⋅2      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqaaewabaaabagaeeysakkaeeoraykaeeiiaaiaeeiiaaiaee4camnaeeyaakmaee4zacmaeeoba4maeiikagiaeeizaqmaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqggpaqkcqggpaqkcqgh9aqpcqqgzbwccqqgpbqacqqgnbwzcqqgubgbcqggoaakcqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapiabcmcapaqaaiabbsfaujabbieaijabbweafjabb6eaojabbccagiabbccagiabbgeabnaabaaaleaacqqgrbwaaeqaaogaeyypa0zaasaaaeaacqggoaakcqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasmaey4kasiaegymaedabeaakiabcmcapiabgucariabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaeiykakiaeyyxictaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqghsislcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabikdayaaaaeaacqqgfbqrcqqgmbatcqqgtbwucqqgfbqrcqqggaaicqqggaaicqqgbbqqdawgaawcbagaee4aasgabeaakiabg2da9maalaaabawaaewaaeaacqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasmaey4kasiaegymaedabeaakiabcmcapiabgwsixpaalaaabagaeeizaqmaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqggpaqkaeaacqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasmaey4kasiaegymaedabeaakiabcmcapiabgucariabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakcaaiabgucariabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaeyyxic9aasaaaeaacqqgkbazcqggoaakcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwacqghrawkcqaixaqmaeqaaogaeiykakiaey4kasiaeeizaqmaeiikagiaeeidaq3aasbaasqaaiabbugarbqabagccqggpaqkaaaacagloagaayzkaagaeyyxictaeiikagiaeeidaq3aasbaasqaaiabbugarjabgucariabigdaxaqabagccqghsislcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabikdayaaaaagaaczcaiaaxmaadaqadaqaaiabisda0agaayjkaiaawmcaaaaa@d0ee@ the gene x is considered differentially expressed in t vs c if the following inequality holds: a >θa     where θa is a threshold to be determined, in correspondence to a significance level α, based on the null hypothesis distribution of a, i.e. the distribution of a derived from experimental replicates. method 3 for each gene x, the expression profiles c and t are fitted using natural cubic splines. the null hypothesis is that both c and t time series share the same fit, the alternative hypothesis is that the fits are not equal. as a first step of the method, the same spline function is fitted simultaneously to the two profiles c and t and the sum of squares residuals ss is calculated; then two different spline functions are fitted separately to time series c and t and the sum of squares residuals ss is computed. to assess differentially expressed genes, the goodness of model fit under the null hypothesis is compared to that under the alternative hypothesis, by calculating a statistic f as: f=ss0−ss1ss0      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqggbgrcqgh9aqpdawcaaqaaiabbofatjabbofatnaacaaaleqabagaegimaadaaogaeyoei0iaee4uamlaee4uam1aawbaasqabeaacqaixaqmaaaakeaacqqgtbwucqqgtbwudaahaawcbeqaaiabicdawaaaaagccawljagaaczcamaabmaabagaegonaydacagloagaayzkaaaaaa@3e08@ gene x is considered differentially expressed if the following inequality holds: f >θf     where θf is a threshold to be determined, in correspondence to a significance level α, based on the null hypothesis distribution of f. null hypothesis distribution model both methods and are based on a threshold derived from the null hypothesis distribution of d and a, respectively, obtained using available replicates. we denote as "available replicates" replicates available for a subset of time samples. let's suppose two experimental replicates are available for a generic time tk and a generic experimental condition. by assuming a log-additive error model as in, the log-expression measurement of each gene x in replicates a and b, can be expressed as: {xa=μ+εaxb=μ+εb      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagabeqaauaabeqaceaaaeaacqqg4baedawgaawcbagaeeyyaegabeaakiabg2da9ggaciab=x7atjabgucariab=v7alnaabaaaleaacqqghbqyaeqaaagcbagaeeieag3aasbaasqaaiabbkgaibqabagccqgh9aqpcqwf8oqbcqghrawkcqwf1oqzdawgaawcbagaeeoyaigabeaaaaaakiaawuhaaiaaxmaacawljawaaewaaeaacqai4aaoaiaawicacaglpaaaaaa@450f@ where μ represents the actual gene expression and εa, εb are two realizations of the error variable ε, monitoring both technical and biological variability. the indices for condition and time tk are omitted here because we refer to pair of replicates available for a generic time sample tk. null hypothesis distribution of variable d to quantify the deviation of expression d under the null hypothesis, equation is applied to available replicates as: dh = xa - xb = εa - εb     different distribution models are used to fit the set of dh values obtained by applying equation to all genes and available replicates. the best model is selected based on the goodness of fit and the parameters precision, and is used as the null hypothesis distribution of d to determine θd to be used in equation to determine the threshold in correspondence to a significance level α, method uses a model to fit the observed statistics rather than using quantiles. the reason for this choice is that the lack of a sufficient number of observations from available replicates renders the determination of appropriate thresholds difficult when low significance levels are chosen, as often the case in microarray studies. if a sufficient number of replicates is available to guarantee a good threshold setting at the desired significance level α, it may be preferable to use quantiles. null hypothesis distribution of variable a at least two replicates for each time sample would be necessary to derive a distribution under the null hypothesis from the data. since we address selection of differentially expressed genes in data-poor condition, i.e. a sufficient number of replicates is not available, a monte carlo procedure is used to derive the null distribution of a. first, dh distribution is derived from d values obtained from available replicates as described above. then, b profiles of length m are sampled from dh under the hypothesis that the error at different time samples is independent and identically distributed. subsequently, b values of ah are calculated from these profiles. finally, different distribution models are used to fit the entire set of ah values and the best model is chosen based on goodness of fit and parameter precision. this model is used as the null hypothesis distribution of a to determine θa to be used in equation as for method a model is fitted to the observed statistics rather than using quantiles to determine the threshold in correspondence to a significance level α. null hypothesis distribution of variable f the null distribution of f is obtained using bootstrap. see for details. intensity dependency of error in affymetrix chips it is well known that dh has an intensity dependent distribution. in particular, analysis of technical replicates of affymetrix human chip has shown that the standardized variable sh: sh0=dh0sddh0=xa−xbsddh0;x¯=xa+xb2      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqabeqacaaabagaee4cam3aawbaasqabeaacqqgibascqaiwaamaagccqgh9aqpdawcaaqaaiabbsgaknaacaaaleqabagaeeisagkaegimaadaaagcbagaee4uamlaeeiraq0aasbaasqaaiabbsgakjabbieaijabbcdawaqabagccqggoaakcuqg4baegaqeaiabcmcapaaacqgh9aqpdawcaaqaaiabbiha4naabaaaleaacqqghbqyaeqaaogaeyoei0iaeeieag3aasbaasqaaiabbkgaibqabaaakeaacqqgtbwucqqgebardawgaawcbagaeeizaqmaeeisagkaegimaadabeaakiabcicaoiqbbiha4zaaragaeiykakcaaiabcuda7aqaaiqbbiha4zaaragaeyypa0zaasaaaeaacqqg4baedawgaawcbagaeeyyaegabeaakiabgucariabbiha4naabaaaleaacqqgibgyaeqaaagcbagaegomaidaaaaacawljagaaczcamaabmaabagaegymaejaegimaadacagloagaayzkaaaaaa@5f0b@ has an intensity independent distribution. therefore, in case of data showing intensity dependency of the variable dh it is convenient to model sh distribution, as indicated in the additional file to derive the threshold θd to be used in equation consistently, the values of d observed from the data are standardized before applying equation if method is used: s=dsddh0);x¯=xt+xc2      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqabeqacaaabagaee4camnaeiikagiaeeidaq3aasbaasqaaiabbugarbqabagccqggpaqkcqgh9aqpdawcaaqaaiabbsgakjabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakcabagaee4uamlaeeiraq0aasbaasqaaiabbsgakjabbieaijabicdawaqabagccqggoaakcuqg4baegaqeaiabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaeiykakcaaiabcuda7aqaaiqbbiha4zaaragaeiikagiaeeidaq3aasbaasqaaiabbugarbqabagccqggpaqkaagaeyypa0zaasaaaeaacqqg4baedaahaawcbeqaaiabbsfaubaakiabcicaoiabbsha0naabaaaleaacqqgrbwaaeqaaogaeiykakiaey4kasiaeeieag3aawbaasqabeaacqqgdbwqaagccqggoaakcqqg0baddawgaawcbagaee4aasgabeaakiabcmcapaqaaiabikdayaaacawljagaaczcamaabmaabagaegymaejaegymaedacagloagaayzkaaaaaa@6586@ analogously, if method is used on data showing intensity dependency of the variable dh ah and θa are derived using sh and the values of a observed from the data are calculated using s instead of d. threshold setting once the null hypothesis distribution of d, a and f are obtained, thresholds θd, θa and θf are determined in correspondence to a significance level α. rather than fixing it a priori, α can be optimized based on a variety of criteria aiming to control the family wise error rate, or the false discovery rate or a compromise between false positive and false negative classification. as an example, let's focus on a criterion based on the control of fdr, defined as the expected proportion of false positive classification among the number sα of genes selected as differentially expressed, using significance level α: fdr=e      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqggbgrcqqgebarcqqgsbgucqgh9aqpcqqgfbqrdawadaqaamaalaaabagaeeoraykaeeiuaafabagaee4uam1aasbaasqaaggaciab=f7ahbqabaaaaagccaglbbgaayzxaagaaczcaiaaxmaadaqadaqaaiabigdaxiabikdayagaayjkaiaawmcaaaaa@3e04@ in case of numerous sets as for microarrays, fdr is well approximated by fdr≈esα      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqqggbgrcqqgebarcqqgsbgucqghijyudawcaaqaaiabbweafjabcufabjabdaeagjabdcfaqjabc2fadbqaaiabdofatnaabaaaleaaiigacqwfxoqyaeqaaaaakiaaxmaacawljawaaewaaeaacqaixaqmcqaizawmaiaawicacaglpaaaaaa@3f45@ calculating fdr requires the estimate of the expected number of false positives, obtained as the product of α by the number n of non differentially expressed genes: e = n · α     n is unknown and is estimated using the bootstrap procedure described in. fdr is calculated for a range of significance levels α and the significance level that guarantees the desired fdr is then used to select differentially expressed genes. since method applies m tests for each gene, the significance level α for method is corrected by applying Šidák correction in order to account for multiple testing. simulation three different experimental conditions with a number of time samples m = were simulated. synthetic data sets were generated for each experimental condition, each consisting of profiles: simulated as differentially expressed and the remaining as random noise. in both cases, the deviation of expression in t vs c was generated at each sampling time tk as standardized deviation: s ~ n ∀k =..., m     where σ was set equal to and μk = for not differentially expressed genes; while, for differentially expressed genes, plausible profiles were obtained by modeling μk as dependent on μk- according to a first order markov model, with the only constraint of being greater than for at least one time samples. samples k = were generated twice for each gene, so as to provide replicates useful to apply method and these replicates were included also in the analysis for method simulated data were used to test the performance of the methods in different experimental conditions. after the null hypothesis distributions of variables d, a and f are modeled as described above, a significance level α had to be fixed to determine the confidence thresholds θd, θa and θf to be used in equations and respectively. we compared the performance of the three methods across the entire range of significance level α by using precision vs recall, and curves of observed false positives divided by the number of selected genes vs number of selected genes. moreover, we compared the average results obtained with the methods by setting α in correspondence to a desired fdr of, in terms of number of selected genes and observed false positives divided by the number of selected genes. all measurements were averaged across the simulations. insulin case study to better appreciate some characteristics of methods and related to the experimental error modeling, the analysis of the null hypothesis distribution of the variables d and a was applied to a real case study on rat muscle cells treated with insulin. the study was performed in vitro, on muscle l rat cell line. cells were treated with insulin at time 0+, just after the collection of a first baseline sample at time 0; eight samples were harvested every hour during eight hours insulin stimulation. a control experiment was also performed in order to be able to distinguish between insulin effect and biological processes of different nature, which take place in the culture simultaneously to insulin induced processes. a total of twenty affymetrix chips rg_u34a were hybridized using four replicates of the basal sample, eight samples collected from the control culture, and eight samples collected from the treated culture. standard affymetrix mas software was used for data pre-processing.  simulation an example of simulated data with m = is shown in figure as clusters obtained by using self organizing maps and pearson correlation as similarity measure on one of the simulated data-set. the average profile is shown for each cluster, together with standard deviation bars; the number of genes for each cluster is also reported. a variety of profiles are represented, such as genes differentially expressed in one or few peaks apparently uncorrelated, profiles that show characteristic bumps and waves of different length or consistent trends along the time series, accompanied by few or numerous time samples with average absolute value greater than the error standard deviation σ. results obtained by applying methods and to simulated data sets are shown as average precision at different ranges of recall intensities in figure left panels. average and standard deviations across the simulations are reported. method slightly outperforms method for precision higher than, since it provides higher values of both precision and recall in all the different sampling conditions. for values of precision lower than, the two curves are almost superimposable. both methods and work better than method for short time series, while for long time series the performance of the methods and is similar, with method slightly superior to methods the areas under the curves are also reported in figure. to better appreciate the ability of the three methods to select differentially expressed genes at low fdr, figure shows the number of observed false positives divided by the number of selected genes for different numbers of selected genes in the range 0–300: method provides number of false positives divided by the number of selected genes always lower than that obtained with method 1; method is slightly superior to method for m = but inferior to both methods and for m = and m = the performance of the three methods with respect to each other is not affected by the number of differentially expressed genes. insulin case study error in this section we present the estimate of the null hypothesis distribution of variable d and variable a. figure shows the values of dh and its standardized value sh plotted vs the average intensity of gene expression values. figure shows the average variance of dh and sh vs the intensity range of gene expression discretized in intervals of constant size. results confirm that while dh has an intensity-dependent distribution, sh has not. the entire set of sh values was fitted using different distribution models. the gaussian mixture model with n = was chosen as the best model. model parameters and their precision are shown in table monte carlo procedure was then used to derive the null distribution of a: b profiles of length m = were sampled from sh distribution and b = values of ah were calculated from these profiles. three different distribution models were used to fit the entire set of ah values. the best fit for ah was obtained with gamma distribution. gene selection we applied the three methods to the data, adopting the false discovery rate as criterion for threshold setting. table shows the number of probe-sets selected by each of the three methods for fdr equal to,,, respectively. results confirm those obtained by using simulation, i.e. the ability of method to select a higher number of differentially expressed genes at controlled fdr for short time series. since focus here is methodology, biological results are not discussed further; confirmation studies and biological interpretation will be presented in a different article. discussion in this work we evaluated the performance of two selection methods here proposed to be applied in time series studies in data-poor conditions, i.e. when the number of available replicates does not make possible or practical the use of standard statistical methods. we also tested the two methods in comparison with a third method from the literature. method compares samples time by time using a statistically based fold change threshold derived from a null hypothesis distribution of variable d. to this purpose, replicates of at least one time sample are necessary. since the threshold is derived based on the experimental variability, method accounts for the error characteristics, e.g. its intensity dependence; therefore, for example using affymetrix chips, genes expressed at high intensities are not penalized with respect to genes at low intensities, which show a higher variability. although method improves upon the use of a constant empirical fold change threshold, it considers time samples independently to each other, which is not a realistic assumption in time series studies. method calculates the area of the region bounded by the time series expression profiles to be compared and considers the gene differentially expressed if this area exceeds a threshold based on a model of the experimental error; therefore, besides accounting for the error distribution, it considers the entire expression profile and not single time samples. also this method needs replicates for at least one time sample in order to derive the null hypothesis distribution. both methods and assumes as working hypothesis that the error at different time points is independent and identically distributed. this hypothesis, on our experience, is usually verified on real data. however, for some experimental settings there may be a dependency of the error on time. in this latter case, it would be more appropriate to perform replicates at different time samples, covering the duration of the experiment, and use method with time dependent threshold settings. methods and can be applied to compare time series from different experimental conditions or a time series vs its baseline, by defining the deviation of expression of gene × for each sample tk as the deviation between treated at time tk and the baseline. moreover, if the sampling grid is different for the two time series, method can be easily generalized by generating ah distribution from b time series with appropriate sampling grid. method uses natural cubic splines to fit time series expression data using a null hypothesis and an alternative hypothesis model, and performs a statistical test on the sum of squares residuals obtained using the two models to assess differential expression. it therefore, besides considering the entire expression profiles, implicitly considers time dependencies. moreover, it does not need any experimental replicate, which makes it practical for a wide range of time series microarray experiments. we tested the performance of the three methods using synthetic data to assess their validity over a range of experimental conditions, specifically the length of the analyzed time series. in the simulation, we used a markov model to generate plausible data, accounting for the dependencies of time samples and not biased toward one of the methods. taking for example μk equal to an arbitrary constant for a random subset of time samples tk in equation would not have accounted for time dependencies, and would have generated non realistic oscillation in the profiles, thus penalizing method which is based on a model fit, with respect to the other methods, in particular method which applies a threshold on each time sample. moreover, simulated profiles represent a variety of possible situations in time series expression, such as profiles characterized by one or few peaks, waves of different length or consistent trends along the time series. results on simulated data showed that method outperforms method independently from the length of the time series being analyzed, probably because, as method is based on single sample comparisons, it is particularly sensitive to random fluctuations due to the noise, thus resulting in a larger number of false positives. method constitutes an improvement with respect to method since the entire expression profile is considered simultaneously and this allows better distinguishing between consistent differences in expression profiles and random oscillations, thus resulting in a lower false negative rate. method as method considers the entire expression profile, but performs better than method only for long time series. looking at the ability of the methods to classify profiles with particular characteristics, we observed that method works better than method to detect differentially expressed genes that show just one or two peaks as differentially expressed. on the opposite, method works better than method in detecting profiles characterized by waves of length greater than three samples, so as profiles that show a characteristic increasing/decreasing trend. method as method is better in detecting bumps and consistent trends in the profiles, than in detecting isolated peaks. however, it needs long time series to perform better than method probably because the fit are more reliable when performed on long time series; method was in fact proposed by the authors on real case studies with more than samples. these results are certainly of interest to address the choice of the algorithm to be used in data-poor time series expression studies, depending on the availability of replicates and on the length of the time series.  microarray time series studies are essential to understand the dynamics of biological molecular events. in order to limit the analysis to those genes that change expression over time, it is necessary to select differentially expressed transcripts. due to the high cost of microarrays, experiments are often performed without replication; therefore, traditional statistical methods can't be applied. here we evaluate the performance of two selection methods applicable in data poor conditions, based on: a statistically based threshold on individual samples; a statistically based threshold to be applied on the area of the region bounded by the time series expression profiles to be compared. application on a real data set on insulin regulation on muscle cells, obtained using affymetrix chips, revealed as the error analysis performed using methods and may be useful to detect error characteristics such as intensity dependencies and to properly address these feature by standardization. we evaluated methods and performance using simulated data with a different number of available samples and compared these performance with those obtained using method based on a splines fit of time series profiles. the results outlines that the two error based methods and work better than method with short time series experiments, while method works better than methods and with long time series experiments. these results might help to optimize the choice of the algorithm to be used in different experimental conditions. a preliminary version of method implemented in r, is available at. authors' contributions bdc conceived the study and performed data analysis under the guidance and supervision of gt. skn was responsible for the overall coordination of microarray experiments realization. ljg performed microarray experiments. cc was responsible for the overall project coordination. all authors read and approved the final manuscript. supplementary material additional file 1 the.pdf document contains a detailed description of the procedure used to standardize the variable dh so as to account for the error intensity dependent distribution, and a detailed description of the procedure used to simulate the data, based on a first order markov model. click here for file acknowledgements this study was supported by ministero dell'università e della ricerca scientificae tecnologica, and by national institute of health, grant roidk this article has been published as part of bmc bioinformatics volume supplement 2007: italian society of bioinformatics: annual meeting the full contents of the supplement are available online at. figures and tables all measurements were averaged across the simulations; standard deviations are reported in parenthesis. number of selected genes and number of false positives divided by the number of selected genes, obtained by applying the methods on simulated data. number of genes selected applying the methods on real data.
method calculates the area of the region bounded by the time series expression profiles, and considers the gene differentially expressed if the area exceeds a threshold based on a model of the experimental error. in order to limit the analysis to those genes that change expression over time, a first necessary step is to select differentially expressed transcripts. application of the three methods to synthetic data indicates that method outperforms the other two both in precision and recall when short time series are analyzed, while method outperforms the other two for long time series. in this data-poor context, we evaluate the performance of three selection methods, using synthetic data, over a range of experimental conditions. method uses a threshold on individual samples based on a model of the experimental error. these two methods are compared to method recently proposed in the literature, which exploits splines fit to compare time series profiles. these results help to address the choice of the algorithm to be used in data-poor time series expression study, depending on the length of the time series. application to real data is also discussed. microarray time series studies are essential to understand the dynamics of molecular events. a variety of methods have been proposed to this purpose; however, these methods are seldom applicable in practice since they require a large number of replicates, often available only for a limited number of samples. three methods are considered, to assess differentially expressed genes in data-poor conditions.
 microarray time series studies are essential to understand the dynamics of molecular events. in order to limit the analysis to those genes that change expression over time, a first necessary step is to select differentially expressed transcripts. a variety of methods have been proposed to this purpose; however, these methods are seldom applicable in practice since they require a large number of replicates, often available only for a limited number of samples. in this data-poor context, we evaluate the performance of three selection methods, using synthetic data, over a range of experimental conditions. application to real data is also discussed.  three methods are considered, to assess differentially expressed genes in data-poor conditions. method uses a threshold on individual samples based on a model of the experimental error. method calculates the area of the region bounded by the time series expression profiles, and considers the gene differentially expressed if the area exceeds a threshold based on a model of the experimental error. these two methods are compared to method recently proposed in the literature, which exploits splines fit to compare time series profiles. application of the three methods to synthetic data indicates that method outperforms the other two both in precision and recall when short time series are analyzed, while method outperforms the other two for long time series.  these results help to address the choice of the algorithm to be used in data-poor time series expression study, depending on the length of the time series. 28– april italian society of bioinformatics: annual meeting bologna, italy 