BACKGROUND
multiplex gene expression studies that are based on microarrays and next-generation sequencing result in the generation of large datasets.
the simultaneous analysis of the expression of thousands of genes in these high-throughput approaches challenges statistical methods and data interpretation.
traditional statistical tests and analysis tools are therefore often insufficient.
one of the most popular types of biological experiments is a two-sample comparison.
gene expression studies often seek to identify genes that are differentially expressed  between rna samples from two types of biological conditions, such as gene knockout mice compared to wild-type mice.
de genes can give insights into biological mechanisms or pathways, and form the basis for further experiments.
the traditional statistical method for identifying de genes between two samples is the student’s t-test.
but a microarray comparative experiment is faced with simultaneously assessing a large number of genes based on a small number of biological or technical replicates.
assessing such a large dataset with a small sample size challenges statistical methods: dealing with many genes but few replicates may lead to large fold changes  driven by outliers, and to small error variances  <cit> .
in order to overcome such problems, the t-test has been modified for microarray data analysis.
the general idea behind these modifications is to obtain more stable estimates of the error variance of a given gene by “borrowing” information from all available genes.
this goal is often obtained by applying  bayesian methods.
examples of these modified tests for microarray data are the sam test  <cit> , the regularized t-test  <cit> , and the b-statistic  <cit> , as reviewed  <cit> .
other statistical approaches for microarray data analysis have introduced linear models  <cit> .
these models allow for more flexibility, for instance when comparing more than two samples or introducing additional sources of variation.
such an example is lin et al.
<cit> , where a complex experimental design and research goal are addressed by setting various contrasts in the design of the linear model.
the bioconductor package limma, developed by smyth  <cit> , applies a gene-wise linear model, and allows for the analysis of complex experiments , as well as more simple replicated experiments with two rna samples.
the package has further developed the ideas of lönsstedt and speed  <cit> , which have been reset within the framework of linear models in order to address the challenges of microarray data analysis.
the statistic  to identify de genes is referred to as the moderated t-statistic, denoted by t*  <cit> .
in the calculation of t*, shrinkage of the estimated error variances towards a pooled estimate is obtained through an empirical bayes approach.
for certain biological problems, it is important to rank genes according to their fc or to impose on a gene to attain a predefined fc threshold before calling it de.
in such situations there remains often a disconnect between the concepts of a statistically significant differential expression  and a biologically meaningful differential expression .
even the modified tests may result in genes with small fcs to be considered statistically significant.
in order to integrate these statistical and biological concepts, a gene can be defined as de when it satisfies a p value and a fc criterion  <cit> .
the advantage of this combined approach is that the t-test or any of the modified tests can be combined with an ad-hoc fc criterion.
the disadvantage of an ad-hoc fc criterion is that it does not take into account error variance and therefore offers no statistical confidence about future results.
in addition, depending on the choice of the fc criterion and significance criterion, various interpretations of the same dataset are possible  <cit> .
the moderated t-statistic  <cit>  has been extended by mccarthy and smyth into a new “test relative to a fc threshold”, abbreviated treat  <cit> .
this method assesses formally whether the true differential expression is greater than a predefined fc criterion.
treat offers greater specificity and reproducibility in identifying de genes, compared to the combined approach of statistical test and ad-hoc fc criterion.
treat has been also added to the limma package.
we have previously applied the nanostring digital platform  <cit>  to study the expression of odorant receptor  genes in mice  <cit> .
in contrast to microarray data, where analog levels of fluorescent intensity are measured, nanostring data represent digital readouts of single molecules in the form of probe counts.
these probes contain unique fluorescent bar codes, and rna abundance of up to  <dig> genes can be analyzed in a single reaction in a single tube.
the nanostring technology thus places itself between qpcr and microarrays in terms of throughput level  <cit> .
we have analyzed with nanostring the expression of half of the or gene repertoire of ~ <dig> genes  <cit> .
here, we have explored statistical tools for our nanostring data, and developed a systematic approach for identifying de genes with respect to a given fc threshold.
we explored the moderated t-statistic , the derivation of which was driven by microarray data , versus the classical t-statistic on comparative nanostring experiments .
we found that t* does not show a protective effect  over t on our nanostring data.
but we also wanted to test whether differential expression is greater than a fc threshold.
therefore we used the analysis relative to a fc threshold together with the classical t-statistic in two approaches.
the first approach is similar to treat as published for t*  <cit> , and we refer to it as ttreat.
then we addressed the arbitrary choice of the fc criterion itself, and developed a two-stage approach, ttreat <dig>  we describe the performance of treat, ttreat, and ttreat <dig> on our nanostring data, both in data simulation experiments and on biological data.
the variability of the fc of a gene is inversely related to the expression level of that gene; lowly expressed genes tend to have a greater error in their measured fc levels  <cit> .
these lowly expressed genes can thus more easily reach a certain fc threshold, and the inverse is true for highly expressed genes.
a non-linear model, the limit fold change model , has been developed to identify de genes based on this relation  <cit> .
we have applied a similar lfc model on our nanostring data.
we did not use the lfc model as a tool for identifying de genes but to set appropriate fc thresholds for genes with various ranges of expression levels, in order to avoid a subjective choice of fc criterion.
the fc thresholds that we thus derived were then used in a subsequent analysis relative to a threshold in order to identify the de genes.
we refer to this setup as the running fc model, and illustrate its use and benefit on biological data.
RESULTS
biological data of odorant receptor gene expression
the main olfactory epithelium is located in the nasal cavity of the mouse, and detects volatile chemicals  in the inhaled air.
the sense of smell must detect chemical stimuli with an immense variety in physicochemical properties.
to accommodate this broad recognition, the mammalian olfactory system has evolved a large repertoire of molecular receptors, odorant receptors.
these receptors are expressed by olfactory sensory neurons in the main olfactory epithelium.
in the mouse, ~ <dig> odorant receptors are encoded by distinct genes, which form the largest gene family in its genome.
it is widely accepted that a mature olfactory sensory neuron expresses only one of the ~ <dig> odorant receptor genes.
the molecular and genetic mechanisms that regulate this one receptor - one neuron rule remain unclear, and are the focus of our research.
we have demonstrated the role of the h element  <cit>  and the p element  <cit>  in the regulation of expression of clusters of odorant receptor genes by genetically engineering mouse strains that lack the h element  <cit>  or the p element  <cit> : these are the Δh and Δp strains and the ΔhxΔp double knockout strain  <cit> .
another mouse strain is the Δolfr7Δ strain  <cit> ; by means of chromosome engineering  <cit> , we excised a  <dig>  megabase region on chromosome  <dig> that contains the olfr <dig> cluster consisting of  <dig> or genes  <cit> .
we have also characterized temporal expression patterns of  <dig> odorant receptor genes in adult and aged mice  <cit> .
here we used datasets  <cit>  from a nanostring analysis of  <dig> or genes comparing knockout versus wild-type mouse strains.
specifically, we used nanostring data obtained from six mutant mice of the ΔhxΔp strain  compared to  <dig> control  mice ; these  <dig> mice are in a mixed genetic background, c57bl/6 j × 129/svev.
another dataset was obtained from six mice of the Δolfr7Δ strain, compared to six control  mice ; these  <dig> mice are in pure genetic background, 129/svev.
with nanostring codeset gorilla, we determined the rna abundance for  <dig> or genes from 1 μg rna of whole olfactory mucosa tissue samples.
each lane of a nanostring cartridge represents a different rna sample and mouse.
thus, there are 6– <dig> biological replicates per biological condition, and no technical replicates.
approaches relative to a fc threshold
our novel method ttreat is similar to treat  <cit> .
it is applied to the regular student’s t-statistic, and requires a predefined fc threshold τ. for the two-stage design, ttreat <dig>  a second threshold θ, with θ>τ, is used in a first “stop and go” stage in which it is decided whether a gene is non-de  or whether it can proceed to the next stage .
then, for all the “go” genes, a ttreat test with fc threshold τ is applied.
in the running fc model, a non-linear model for fc versus average gene expression is first used to determine various fc thresholds τi for a number of ranges of expression levels.
genes are then binned in k gene expression levels, and the appropriate τi is used per concentration bin in a subsequent analysis relative to a fc threshold.
simulated data
because nanostring data are not yet as widely studied as microarray data, we have conducted a data simulation procedure that represents one of our typical two-group comparison nanostring experiments.
the procedure for simulating one nanostring dataset was conducted according to the following steps and distributional assumptions:
• to get a general idea about the variances of nanostring gene expression data, the genes in the ΔhxΔp dataset were used as an example gene population.
biological data from ΔhxΔp mice were chosen, as they represent a noisier dataset due to the mixed genetic background of this strain  <cit> : the resulting simulated data will not represent the cleanest example.
the ΔhxΔp dataset was used only for the next step of the simulation exercise.
• subsequently,  <dig> real variances across genes, σg <dig>  were drawn from an inverse-gamma distribution: inv − gamma.
the parameters g <dig> and g <dig> are estimated by a maximum likelihood procedure in the above described gene population.
• each of the  <dig> σg <dig> gave rise to three randomly drawn real differences βg.
the σg <dig> and three corresponding βg were included in one of the following three gene groups:
de genes: the βg were drawn from a gaussian ~n <dig> σgvstart and had to satisfy the criterion: |βg| ≥ log <dig>
non-de genes, noisy: similarly, βg  drawn from ~n <dig> σgvstart but with |βg| < log <dig>
non-de genes, regular: obviously for this group the real βg were set to  <dig>
note that empirically, the normal distribution seems acceptable for the βg, based on our nanostring data.
• the βg that were produced in the previous step served as true differences that were then subsequently used to simulate three possible estimates β^g from a gaussian ~nβg,σgzg.
similarly residual variances were drawn from a chi-square distribution:
~σg2dfgχ2dfg
since the quantity ω used to initiate the simulation as described above defines the de genes by the rule |βg| ≥ log <dig>  the distributional center of the β^g of de genes actually lies a little further than the actual value of ω, depending on the value of their variance: σgzg.
henceforth, we will refer to ω as the fc ω with respect to which the data was simulated.
every such simulated dataset was thus initially based on  <dig> σg <dig> leading to a total of  <dig> genes.
the regular non-de genes  were fixed as 20% of the  <dig> genes.
however, to assess whether the number of de genes in a dataset influences the results, the percentage of de genes  was varied between 1% and 40% .
on the simulated data, the following statistical approaches for identifying genes that are de with regard to a certain fc threshold were assessed:
treat with regard to fc threshold τ
ttreat with regard to fc threshold τ
ttreat <dig> with a bilateral p value calculation in the second stage with regard to fc thresholds θ  and τ
as the two-stage design ttreat <dig> relies on the choice of two related fc thresholds, one for each of the two stages, it is not straightforward to compare the performance of ttreat <dig> to treat and/or ttreat.
therefore we have derived various experimental schemes to show in which situations the use of a two-stage design like ttreat <dig> can be beneficial.
results on simulated data
we simulated  <dig> realizations of a two-group comparison nanostring experiment with  <dig> genes for which the de genes are simulated with respect to a certain fc difference ω as described above.
as such we know beforehand to which group  each gene in these  <dig> datasets belongs.
we can therefore assess the performance of the procedures described in the methods section.
for this purpose, we fix vstart =  <dig> and the significance  was set at  <dig>  or at  <dig> .
the mean of the false discoveries  and missed genes  over the  <dig> generated datasets is plotted for the various statistical approaches for fcs ω  and τ .
we also report mean area under the roc-curve  values over the  <dig> generated data sets.
figure  <dig> and additional file  <dig> compare treat and ttreat on a simulation exercise where the fc difference ω and the fc threshold τ were chosen to be equal.
we tested values for ω = τ of  <dig>  ,  <dig>   and  <dig> .
across the various percentages of de genes , the percentages next to the gray and black bars of ttreat in figure  <dig> represent the percentual decrease  or increase  in false discoveries or missed genes with respect to the reference, treat.
we find that ttreat always results in fewer false discoveries than treat, at p =  <dig>  and p =  <dig> .
moreover, when p =  <dig>  is applied as a significance cut-off , ttreat decreases the false discoveries for only a few more missed genes compared to treat .
for p =  <dig>   the latter is true for datasets with small percentages of de genes but for datasets with ~10% of de genes, the decrease in false discoveries is nullified by the increase in missed genes .
in terms of auc performance there is no difference between treat and ttreat across the percentages of de genes .
simulated fc between test and control =  <dig>
simulated fc between test and control =  <dig>
simulated fc between test and control = 2
simulated fc between tets and control =  <dig>  and stringent test
simulated fc between test and control =  <dig>  and non-stringent test
*mean auc over the  <dig> simulated data sets with 95% prediction intervals.
in order to illustrate the added value of the safety margin around the fc threshold τ that has been chosen for statistical analysis, two simulation experiments were done.
in the first experiment, the data were simulated such that the actual de genes had a fc difference ω of at least  <dig>  , and p was set at  <dig> .
a blind, stringent ttreat test at the higher fc threshold τ of  <dig>   will result in a high number of missed genes .
the reference test for this scenario is a ttreat with fc threshold τ = ω =  <dig> , represented by the black and gray stacked bars in figure 2b.
.
the stringent ttreat  decreases, as expected, the false discoveries, by 90% in comparison to the reference test  but at the cost of a 50% increase in missed genes .
when a two-stage approach is used with fcs θ and τ chosen at  <dig>  and  <dig> , the high number of missed genes of the blind, stringent ttreat test  is reduced .
at the same time, the two-stage approach has a protective effect over the reference ttreat test  by decreasing the number of false positives .
the stringent test has a much lower auc than the reference test and two-stage approach .
the two-stage approach has the best overall performance  for this experiment.
a second simulation experiment focuses on the inverse situation.
data are simulated with respect to a fc of ω =  <dig>  , and p =  <dig> .
one ttreat test aims too low with a fc threshold of  <dig>  , leading to many false discoveries.
the reference test for this scenario is a ttreat with fc threshold τ = ω =  <dig> , represented by the black and gray stacked bars in figure 3b.
.
a two-stage approach decreases the number of false positives that are obtained by the non-stringent test with too low a fc threshold  while maintaining the false negatives .
the auc values show that the overall performance of these three tests is very similar .
benefit of ttreat <dig> for a biological dataset
we analyzed samples of six Δolfr7Δ mutant mice in comparison to six control mice with our nanostring gorilla codeset containing  <dig> or genes.
the olfr <dig> cluster contains  <dig> or genes, and for  <dig> or genes we could design specific probes.
of these  <dig> genes,  <dig> were defined as informative, as the median of the normalized nanostring counts in the control mice is at least  <dig>  because these  <dig> or genes are not present in the genome of the mutant mice, the counts for these or genes in mutant mice should be at background levels, and these genes are a priori de.
the ∆olfr7∆ strain is thus the ultimate negative control for probe specificity.
because the maximum counts for the negative controls is ~ <dig>  we reasoned that even for expressed olfr <dig> cluster genes with low normalized counts , we could expect a >2-fold change.
we thus set the fc threshold τ at  <dig> in a first ttreat analysis, with p =  <dig> .
the mc plot in figure 4a shows the results: of the  <dig> genes in the deleted olfr <dig> cluster,  <dig> are detected as having a fc significantly lower than ½, but one deleted gene, olfr <dig>  is missed.
we lowered the threshold to τ =  <dig>  in a second ttreat analysis: a total of  <dig> genes is now identified with a fc significantly lower than 1/ <dig>  or higher than  <dig>  .
but two genes with a fc higher than  <dig>  turn out to reside outside the deleted cluster: olfr <dig>  which could perhaps make sense given it resides close to the deleted region, and olfr <dig>  which is on chromosome  <dig> and ought not to be affected.
importantly, when we apply a two-stage approach, ttreat <dig> with θ =  <dig> and τ =  <dig> , the  <dig> genes of the deleted cluster, and only these genes, are found to be de : no genes are missed, and no genes outside the cluster are identified as de.
the experimentally rare case of a deletion mutant, in which certain genes are not present, thus enables us to demonstrate the benefit of ttreat <dig>
the running fc model
we applied the running fc model on our nanostring data obtained from six ΔhxΔp mutant mice and  <dig> control mice.
as these mice are from a mixed genetic background, we used a less-stringent ttreat with τ =  <dig> , and p =  <dig> .
the results are shown as an ma plot in figure 5a and as an mc plot in figure 5c.
we find that nine genes from the p cluster and  <dig> genes from the h cluster have a fc significantly lower than 1/ <dig> .
but two additional genes, olfr <dig> on chromosome  <dig> and olfr <dig> on chromosome  <dig>  are identified as de by ttreat with τ =  <dig> .
when we apply a running fc model, a fc threshold value τ is not chosen.
instead the model calculates a different τ for ten gene expression levels.
these values are ≥ <dig>  for the lower expression levels and drop to  <dig>  for the very high expression levels .
after calculating the τ values, we applied them together with a ttreat test in our running fc model.
this model identified  <dig> de genes,  <dig> of which were also identified by ttreat with τ =  <dig>  .
the lowly expressed gene olfr <dig>  is no longer identified as de by the running fc model, but the model identified an additional or gene, olfr <dig>  which is plausible as it resides within the p cluster.the advantage of tests relative to a fc threshold over combined approaches  is illustrated in figure 5a.
by combining the regular t-test and a fc criterion of  <dig> , a total of  <dig> genes are identified as de.
of these,  <dig> reside outside the h and p clusters , and are thus in reality most likely not de.
discussion
gene expression has been studied extensively for a wide variety of biological problems.
over the years various techniques have been developed for and applied to gene expression analysis.
ordered from low to high throughput, these techniques include qpcr, nanostring, microarrays, and rna sequencing.
these techniques come with their own specifications, advantages and disadvantages  <cit>  but each also with specific challenges for statistical analysis.
microarrays have thus far been the most widely used approach to study gene expression, and have therefore led to numerous analytical tools and statistical tests.
many of these tools can also be applied to nanostring data or qpcr data, and some have resulted in practical software packages for nanostring  <cit> .
the moderated t  based treat, which was developed for microarray data  <cit>  works well on our nanostring data.
but we found that ttreat with regular t is even more appropriate.
moreover, we propose an alternative approach that defines a safety margin around the fc threshold.
this ttreat <dig> application in two subsequent stages can be beneficial when data are expected to be noisy around a chosen fc threshold.
we also developed a technique for finding fc thresholds that vary with the expression level of genes, followed by applying these objectively set thresholds in ttreat.
in principle, this running fc model can be applied to gene expression data obtained with any technique.
the goal of this study is not to determine which test has the highest performance on nanostring data, although we describe the overall performance  of treat, ttreat and ttreat <dig> on simulated data.
statistical tests like treat, fine-tuned to gene expression data, are difficult to outperform.
in certain situations, however, an alternative test may prove more beneficial for the technology that is used and the problem that is studied.
in our studies of odorant receptor gene expression, a falsely discovered gene leads to extensive follow-up experiments such as in situ hybridization and even generation of a new mouse strain by gene targeting.
such follow-up experiments are costly, lengthy, and time-consuming.
we thus need to minimize false discoveries.
when using nanostring data, we find that ttreat protects more against false discoveries at significance cut-offs of p =  <dig>  and p =  <dig> .
it is perhaps counterintuitive that ttreat instead of treat is more beneficial for our dataset, as treat and the moderated t statistic are known for their ability to decrease false discoveries in microarray experiments.
there are four main messages.
first, we demonstrate with data simulations that ttreat is more appropriate for our nanostring data than treat.
we used biological data to initiate data simulation experiments, and find that ttreat results in fewer false discoveries compared to treat when using significance cut-offs of p =  <dig>  or p =  <dig> .
but when the percentage of de genes in the simulated data is high, the advantage of ttreat over treat in terms of false discoveries is countered by ttreat missing more de genes.
when the overall performance of treat and ttreat is compared by aucs , there is no difference on our simulated data.
the treat bayesian approach shrinks  the individual gene sample variances towards a pooled estimate, which works particularly well when the number of replicates  is small; often no more than two replicates are available for microarrays  <cit> .
the nanostring counts in our dataset are not noisy and the throughput is much lower  compared to a microarray experiment.
so some of the gene-wise variances may be shrunk  by too large a factor when applying treat on our nanostring data, as the hierarchical bayesian model is expecting a few outlier variances.
second, when noisy data are expected, with many non-de genes showing signs of being increased or decreased, it is beneficial to use a safety margin around the chosen fc threshold in an analysis relative to a fc threshold.
a good choice seems to set θ to τ +  <dig>  or τ +  <dig>
third, the particular choice of a fc threshold can influence the interpretation of biological data  <cit> .
the running fc model can be applied in order to avoid the subjective choice of a fc threshold, or when data vary strongly with expression levels.
we used the same dataset to determine the fc thresholds for various gene expression levels and to apply ttreat, which may seem suboptimal.
it would be ideal to determine fc thresholds on an independent dataset before using these thresholds in the subsequent analysis of the actual data.
but often obtaining these independent data may be too expensive or impossible.
we believe that performing the running fc model on the same dataset gives accurate results despite its data-driven nature in such circumstances.
fourth, the choice of analytical tools must be driven by the research goals, the gene expression technique, and the hypothesis that is being tested.
treat, ttreat, or ttreat <dig> either give a high number of false discoveries and few missed genes, only a few false discoveries for a higher number of missed genes, or an average but similar number of false discoveries and missed genes.
the inevitable trade-off between false discoveries and missed genes must be evaluated dependent on the experimental context.
in our case - nanostring studies of odorant receptor gene expression - we want to minimize the false discoveries, because the validation and follow-up experiments are laborious and expensive.
but for other projects, the inverse may be desirable: the number of missed genes needs to be minimized, for instance in large-scale screening of molecules against a drug target, where missing a potential blockbuster drug cannot be afforded.
in the running fc model, the choice of the fc percentile in each bin in the development of the lfc model  as well as the choice of the k bins when applying the model  affect the final results of the running fc model, thus the de genes discovered.
depending on the percentage of de genes that is expected, percentile values ranging from p <dig>  to p <dig>  are deemed acceptable.
regarding the number of k bins, a good coverage of the full expression range of the genes under analysis should be ensured.
performance results on a single simulated dataset  would have been subject to chance .
repeating the simulations  <dig> times permits an understanding of the variability in test performance that may be obtained between different simulations.
the variability of our mean results is assessed by the 95% prediction intervals around the mean number of false discoveries  on these  <dig> data simulations.
neither of the statistical tests shows outliers in performance that could contradict our conclusions.
our data simulation procedure requires a “real” data input at the beginning only to get an idea about the distribution of gene error variances for nanostring.
this input came from our data as described in the results.
we also used an independent, published nanostring data set with ~ <dig> mouse genes   <cit>  to initiate the set of data simulation experiments again.
we find that simulations initiated by these data lead to very similar results .
we present approaches for testing the significance of gene expression relative to a fc threshold, with a focus on our nanostring data on or genes.
when carrying out these hundreds of significance tests simultaneously, the issue of multiple testing can be raised.
a significance analysis for microarrays  test has been developed  <cit>  that encompasses both a test statistic for ranking  and a way to control the false discovery rate and thus to adjust for multiple testing.
the advantage of our approaches is that they provide p values that can be adjusted by a multiple testing method of choice, ranging from bonferroni to benjamini and hochberg  <cit> .
there is a wide variety of adjustment procedures for multiple testing.
the choice of the most appropriate procedure for a given analysis is not straightforward.
guidelines and criteria to be attributed to the multiple testing approaches have been described  <cit> .
in general, the formulation of a composite null-hypothesis forces treat, ttreat, and ttreat <dig> to rely on a more conservative way of p value calculation  <cit> .
the distribution of the resulting p values is thus skewed towards the larger values, which has consequences:  if a method for adjusting multiple testing problems relies on the uniformity of the distribution of the p values, it cannot be used as this criterion is not met and  we feel that when applied to nanostring data , these already conservative p values should not be adjusted by yet another conservative method.
for the development of treat and initially also for t*  <cit> , it was decided to use gene-wise linear models with arbitrary coefficients and contrasts of interest  as a contextual and practical background; see also the limma package  <cit> .
the ttreat, ttreat <dig>  and the running fc model are also based on gene-wise linear models.
others  <cit>  have used linear models, or, more precisely, analysis of variance models  to assess de genes in microarray experiments.
the main difference is that all genes are included in a single linear model, presenting the advantage that the normalization process is combined with the data analysis and thus pre-normalization of the data is not necessary.
obviously for such an anova model to be justified, its underlying model assumptions should be met, but they can be assessed quite straightforward.
for a smaller nanostring codeset century  <cit> , which consists of  <dig> or genes and  <dig> reference genes, we developed anova models that also dealt with assessing the de genes relative to a fc threshold.
the negative controls could not be included in these anova models, as their residuals were making the general residual distribution heavily tailed and therefore non-gaussian.
equality of variances was not met but within an acceptable borderline range on these genes.
the results of these anova models were not satisfactory, as the genes that were identified as de genes  were sometimes very different from the treat, ttreat and ttreat <dig> results, and less biologically meaningful.
anova models including various genes in the same linear model can be applied on nanostring data and may be very useful for certain purposes.
on our data, however, the gene-wise linear models and thus treat, ttreat, or ttreat <dig> were preferred over an approach that models all genes simultaneously.
CONCLUSIONS
our gene-wise statistical analysis of gene expression data with significance relative to a fc threshold gives reproducible and reliable results on nanostring data of odorant receptor genes, the largest gene family in the mouse genome.
because it is difficult to set a biologically meaningful fc threshold in advance, we have developed methods that provide guidance in determining a stable fc threshold  or in avoiding this choice altogether.
by applying a two-stage approach, a safety margin around the fc threshold can be set, which is beneficial in certain situations.
our running fc model identifies fc thresholds in an automated way, and is thus a more objective model leading to results that are more reproducible.
this model is well suited for the problem of higher fc variances of lower expression values, thus avoiding a bias.

 multiplex gene expression studies that are based on microarrays and next-generation sequencing result in the generation of large datasets. the simultaneous analysis of the expression of thousands of genes in these high-throughput approaches challenges statistical methods and data interpretation. traditional statistical tests and analysis tools are therefore often insufficient. one of the most popular types of biological experiments is a two-sample comparison. gene expression studies often seek to identify genes that are differentially expressed between rna samples from two types of biological conditions, such as gene knockout mice compared to wild-type mice. de genes can give insights into biological mechanisms or pathways, and form the basis for further experiments. the traditional statistical method for identifying de genes between two samples is the student’s t-test. but a microarray comparative experiment is faced with simultaneously assessing a large number of genes based on a small number of biological or technical replicates. assessing such a large dataset with a small sample size challenges statistical methods: dealing with many genes but few replicates may lead to large fold changes driven by outliers, and to small error variances. in order to overcome such problems, the t-test has been modified for microarray data analysis. the general idea behind these modifications is to obtain more stable estimates of the error variance of a given gene by “borrowing” information from all available genes. this goal is often obtained by applying bayesian methods. examples of these modified tests for microarray data are the sam test, the regularized t-test, and the b-statistic, as reviewed. other statistical approaches for microarray data analysis have introduced linear models. these models allow for more flexibility, for instance when comparing more than two samples or introducing additional sources of variation. such an example is lin et al., where a complex experimental design and research goal are addressed by setting various contrasts in the design of the linear model. the bioconductor package limma, developed by smyth, applies a gene-wise linear model, and allows for the analysis of complex experiments, as well as more simple replicated experiments with two rna samples. the package has further developed the ideas of lönsstedt and speed, which have been reset within the framework of linear models in order to address the challenges of microarray data analysis. the statistic to identify de genes is referred to as the moderated t-statistic, denoted by t*. in the calculation of t*, shrinkage of the estimated error variances towards a pooled estimate is obtained through an empirical bayes approach. for certain biological problems, it is important to rank genes according to their fc or to impose on a gene to attain a predefined fc threshold before calling it de. in such situations there remains often a disconnect between the concepts of a statistically significant differential expression and a biologically meaningful differential expression. even the modified tests may result in genes with small fcs to be considered statistically significant. in order to integrate these statistical and biological concepts, a gene can be defined as de when it satisfies a p value and a fc criterion. the advantage of this combined approach is that the t-test or any of the modified tests can be combined with an ad-hoc fc criterion. the disadvantage of an ad-hoc fc criterion is that it does not take into account error variance and therefore offers no statistical confidence about future results. in addition, depending on the choice of the fc criterion and significance criterion, various interpretations of the same dataset are possible. the moderated t-statistic has been extended by mccarthy and smyth into a new “test relative to a fc threshold”, abbreviated treat. this method assesses formally whether the true differential expression is greater than a predefined fc criterion. treat offers greater specificity and reproducibility in identifying de genes, compared to the combined approach of statistical test and ad-hoc fc criterion. treat has been also added to the limma package. we have previously applied the nanostring digital platform to study the expression of odorant receptor genes in mice. in contrast to microarray data, where analog levels of fluorescent intensity are measured, nanostring data represent digital readouts of single molecules in the form of probe counts. these probes contain unique fluorescent bar codes, and rna abundance of up to genes can be analyzed in a single reaction in a single tube. the nanostring technology thus places itself between qpcr and microarrays in terms of throughput level. we have analyzed with nanostring the expression of half of the or gene repertoire of ~ genes. here, we have explored statistical tools for our nanostring data, and developed a systematic approach for identifying de genes with respect to a given fc threshold. we explored the moderated t-statistic, the derivation of which was driven by microarray data, versus the classical t-statistic on comparative nanostring experiments. we found that t* does not show a protective effect over t on our nanostring data. but we also wanted to test whether differential expression is greater than a fc threshold. therefore we used the analysis relative to a fc threshold together with the classical t-statistic in two approaches. the first approach is similar to treat as published for t*, and we refer to it as ttreat. then we addressed the arbitrary choice of the fc criterion itself, and developed a two-stage approach, ttreat we describe the performance of treat, ttreat, and ttreat on our nanostring data, both in data simulation experiments and on biological data. the variability of the fc of a gene is inversely related to the expression level of that gene; lowly expressed genes tend to have a greater error in their measured fc levels. these lowly expressed genes can thus more easily reach a certain fc threshold, and the inverse is true for highly expressed genes. a non-linear model, the limit fold change model, has been developed to identify de genes based on this relation. we have applied a similar lfc model on our nanostring data. we did not use the lfc model as a tool for identifying de genes but to set appropriate fc thresholds for genes with various ranges of expression levels, in order to avoid a subjective choice of fc criterion. the fc thresholds that we thus derived were then used in a subsequent analysis relative to a threshold in order to identify the de genes. we refer to this setup as the running fc model, and illustrate its use and benefit on biological data.  biological data of odorant receptor gene expression the main olfactory epithelium is located in the nasal cavity of the mouse, and detects volatile chemicals in the inhaled air. the sense of smell must detect chemical stimuli with an immense variety in physicochemical properties. to accommodate this broad recognition, the mammalian olfactory system has evolved a large repertoire of molecular receptors, odorant receptors. these receptors are expressed by olfactory sensory neurons in the main olfactory epithelium. in the mouse, ~ odorant receptors are encoded by distinct genes, which form the largest gene family in its genome. it is widely accepted that a mature olfactory sensory neuron expresses only one of the ~ odorant receptor genes. the molecular and genetic mechanisms that regulate this one receptor - one neuron rule remain unclear, and are the focus of our research. we have demonstrated the role of the h element and the p element in the regulation of expression of clusters of odorant receptor genes by genetically engineering mouse strains that lack the h element or the p element: these are the Δh and Δp strains and the ΔhxΔp double knockout strain. another mouse strain is the Δolfr7Δ strain; by means of chromosome engineering, we excised a megabase region on chromosome that contains the olfr cluster consisting of or genes. we have also characterized temporal expression patterns of odorant receptor genes in adult and aged mice. here we used datasets from a nanostring analysis of or genes comparing knockout versus wild-type mouse strains. specifically, we used nanostring data obtained from six mutant mice of the ΔhxΔp strain compared to control mice; these mice are in a mixed genetic background, c57bl/6 j × 129/svev. another dataset was obtained from six mice of the Δolfr7Δ strain, compared to six control mice; these mice are in pure genetic background, 129/svev. with nanostring codeset gorilla, we determined the rna abundance for or genes from 1 μg rna of whole olfactory mucosa tissue samples. each lane of a nanostring cartridge represents a different rna sample and mouse. thus, there are 6– biological replicates per biological condition, and no technical replicates. approaches relative to a fc threshold our novel method ttreat is similar to treat. it is applied to the regular student’s t-statistic, and requires a predefined fc threshold τ. for the two-stage design, ttreat a second threshold θ, with θ>τ, is used in a first “stop and go” stage in which it is decided whether a gene is non-de or whether it can proceed to the next stage. then, for all the “go” genes, a ttreat test with fc threshold τ is applied. in the running fc model, a non-linear model for fc versus average gene expression is first used to determine various fc thresholds τi for a number of ranges of expression levels. genes are then binned in k gene expression levels, and the appropriate τi is used per concentration bin in a subsequent analysis relative to a fc threshold. simulated data because nanostring data are not yet as widely studied as microarray data, we have conducted a data simulation procedure that represents one of our typical two-group comparison nanostring experiments. the procedure for simulating one nanostring dataset was conducted according to the following steps and distributional assumptions: • to get a general idea about the variances of nanostring gene expression data, the genes in the ΔhxΔp dataset were used as an example gene population. biological data from ΔhxΔp mice were chosen, as they represent a noisier dataset due to the mixed genetic background of this strain: the resulting simulated data will not represent the cleanest example. the ΔhxΔp dataset was used only for the next step of the simulation exercise. • subsequently, real variances across genes, σg were drawn from an inverse-gamma distribution: inv − gamma. the parameters g and g are estimated by a maximum likelihood procedure in the above described gene population. • each of the σg gave rise to three randomly drawn real differences βg. the σg and three corresponding βg were included in one of the following three gene groups: de genes: the βg were drawn from a gaussian ~n σgvstart and had to satisfy the criterion: |βg| ≥ log non-de genes, noisy: similarly, βg  drawn from ~n σgvstart but with |βg| < log non-de genes, regular: obviously for this group the real βg were set to note that empirically, the normal distribution seems acceptable for the βg, based on our nanostring data. • the βg that were produced in the previous step served as true differences that were then subsequently used to simulate three possible estimates β^g from a gaussian ~nβg,σgzg. similarly residual variances were drawn from a chi-square distribution: ~σg2dfgχ2dfg since the quantity ω used to initiate the simulation as described above defines the de genes by the rule |βg| ≥ log the distributional center of the β^g of de genes actually lies a little further than the actual value of ω, depending on the value of their variance: σgzg. henceforth, we will refer to ω as the fc ω with respect to which the data was simulated. every such simulated dataset was thus initially based on σg leading to a total of genes. the regular non-de genes were fixed as 20% of the genes. however, to assess whether the number of de genes in a dataset influences the results, the percentage of de genes was varied between 1% and 40%. on the simulated data, the following statistical approaches for identifying genes that are de with regard to a certain fc threshold were assessed: treat with regard to fc threshold τ ttreat with regard to fc threshold τ ttreat with a bilateral p value calculation in the second stage with regard to fc thresholds θ and τ as the two-stage design ttreat relies on the choice of two related fc thresholds, one for each of the two stages, it is not straightforward to compare the performance of ttreat to treat and/or ttreat. therefore we have derived various experimental schemes to show in which situations the use of a two-stage design like ttreat can be beneficial. results on simulated data we simulated realizations of a two-group comparison nanostring experiment with genes for which the de genes are simulated with respect to a certain fc difference ω as described above. as such we know beforehand to which group each gene in these datasets belongs. we can therefore assess the performance of the procedures described in the methods section. for this purpose, we fix vstart = and the significance was set at or at. the mean of the false discoveries and missed genes over the generated datasets is plotted for the various statistical approaches for fcs ω and τ. we also report mean area under the roc-curve values over the generated data sets. figure  and additional file compare treat and ttreat on a simulation exercise where the fc difference ω and the fc threshold τ were chosen to be equal. we tested values for ω = τ of, and. across the various percentages of de genes, the percentages next to the gray and black bars of ttreat in figure  represent the percentual decrease or increase in false discoveries or missed genes with respect to the reference, treat. we find that ttreat always results in fewer false discoveries than treat, at p =  and p = . moreover, when p =  is applied as a significance cut-off, ttreat decreases the false discoveries for only a few more missed genes compared to treat. for p =  the latter is true for datasets with small percentages of de genes but for datasets with ~10% of de genes, the decrease in false discoveries is nullified by the increase in missed genes. in terms of auc performance there is no difference between treat and ttreat across the percentages of de genes. simulated fc between test and control =  simulated fc between test and control =  simulated fc between test and control = 2 simulated fc between tets and control =  and stringent test simulated fc between test and control =  and non-stringent test *mean auc over the simulated data sets with 95% prediction intervals. in order to illustrate the added value of the safety margin around the fc threshold τ that has been chosen for statistical analysis, two simulation experiments were done. in the first experiment, the data were simulated such that the actual de genes had a fc difference ω of at least, and p was set at. a blind, stringent ttreat test at the higher fc threshold τ of will result in a high number of missed genes. the reference test for this scenario is a ttreat with fc threshold τ = ω = , represented by the black and gray stacked bars in figure 2b.. the stringent ttreat decreases, as expected, the false discoveries, by 90% in comparison to the reference test but at the cost of a 50% increase in missed genes. when a two-stage approach is used with fcs θ and τ chosen at and, the high number of missed genes of the blind, stringent ttreat test is reduced. at the same time, the two-stage approach has a protective effect over the reference ttreat test by decreasing the number of false positives. the stringent test has a much lower auc than the reference test and two-stage approach. the two-stage approach has the best overall performance for this experiment. a second simulation experiment focuses on the inverse situation. data are simulated with respect to a fc of ω =, and p = . one ttreat test aims too low with a fc threshold of, leading to many false discoveries. the reference test for this scenario is a ttreat with fc threshold τ = ω = , represented by the black and gray stacked bars in figure 3b.. a two-stage approach decreases the number of false positives that are obtained by the non-stringent test with too low a fc threshold while maintaining the false negatives. the auc values show that the overall performance of these three tests is very similar. benefit of ttreat for a biological dataset we analyzed samples of six Δolfr7Δ mutant mice in comparison to six control mice with our nanostring gorilla codeset containing or genes. the olfr cluster contains or genes, and for or genes we could design specific probes. of these genes, were defined as informative, as the median of the normalized nanostring counts in the control mice is at least because these or genes are not present in the genome of the mutant mice, the counts for these or genes in mutant mice should be at background levels, and these genes are a priori de. the ∆olfr7∆ strain is thus the ultimate negative control for probe specificity. because the maximum counts for the negative controls is ~ we reasoned that even for expressed olfr cluster genes with low normalized counts, we could expect a >2-fold change. we thus set the fc threshold τ at in a first ttreat analysis, with p = . the mc plot in figure 4a shows the results: of the genes in the deleted olfr cluster, are detected as having a fc significantly lower than ½, but one deleted gene, olfr is missed. we lowered the threshold to τ = in a second ttreat analysis: a total of genes is now identified with a fc significantly lower than 1/ or higher than. but two genes with a fc higher than turn out to reside outside the deleted cluster: olfr which could perhaps make sense given it resides close to the deleted region, and olfr which is on chromosome and ought not to be affected. importantly, when we apply a two-stage approach, ttreat with θ = and τ =, the genes of the deleted cluster, and only these genes, are found to be de: no genes are missed, and no genes outside the cluster are identified as de. the experimentally rare case of a deletion mutant, in which certain genes are not present, thus enables us to demonstrate the benefit of ttreat the running fc model we applied the running fc model on our nanostring data obtained from six ΔhxΔp mutant mice and control mice. as these mice are from a mixed genetic background, we used a less-stringent ttreat with τ =, and p = . the results are shown as an ma plot in figure 5a and as an mc plot in figure 5c. we find that nine genes from the p cluster and genes from the h cluster have a fc significantly lower than 1/. but two additional genes, olfr on chromosome and olfr on chromosome are identified as de by ttreat with τ =. when we apply a running fc model, a fc threshold value τ is not chosen. instead the model calculates a different τ for ten gene expression levels. these values are ≥ for the lower expression levels and drop to for the very high expression levels. after calculating the τ values, we applied them together with a ttreat test in our running fc model. this model identified de genes, of which were also identified by ttreat with τ =. the lowly expressed gene olfr is no longer identified as de by the running fc model, but the model identified an additional or gene, olfr which is plausible as it resides within the p cluster.the advantage of tests relative to a fc threshold over combined approaches is illustrated in figure 5a. by combining the regular t-test and a fc criterion of, a total of genes are identified as de. of these, reside outside the h and p clusters, and are thus in reality most likely not de. discussion gene expression has been studied extensively for a wide variety of biological problems. over the years various techniques have been developed for and applied to gene expression analysis. ordered from low to high throughput, these techniques include qpcr, nanostring, microarrays, and rna sequencing. these techniques come with their own specifications, advantages and disadvantages but each also with specific challenges for statistical analysis. microarrays have thus far been the most widely used approach to study gene expression, and have therefore led to numerous analytical tools and statistical tests. many of these tools can also be applied to nanostring data or qpcr data, and some have resulted in practical software packages for nanostring. the moderated t based treat, which was developed for microarray data works well on our nanostring data. but we found that ttreat with regular t is even more appropriate. moreover, we propose an alternative approach that defines a safety margin around the fc threshold. this ttreat application in two subsequent stages can be beneficial when data are expected to be noisy around a chosen fc threshold. we also developed a technique for finding fc thresholds that vary with the expression level of genes, followed by applying these objectively set thresholds in ttreat. in principle, this running fc model can be applied to gene expression data obtained with any technique. the goal of this study is not to determine which test has the highest performance on nanostring data, although we describe the overall performance of treat, ttreat and ttreat on simulated data. statistical tests like treat, fine-tuned to gene expression data, are difficult to outperform. in certain situations, however, an alternative test may prove more beneficial for the technology that is used and the problem that is studied. in our studies of odorant receptor gene expression, a falsely discovered gene leads to extensive follow-up experiments such as in situ hybridization and even generation of a new mouse strain by gene targeting. such follow-up experiments are costly, lengthy, and time-consuming. we thus need to minimize false discoveries. when using nanostring data, we find that ttreat protects more against false discoveries at significance cut-offs of p =  and p = . it is perhaps counterintuitive that ttreat instead of treat is more beneficial for our dataset, as treat and the moderated t statistic are known for their ability to decrease false discoveries in microarray experiments. there are four main messages. first, we demonstrate with data simulations that ttreat is more appropriate for our nanostring data than treat. we used biological data to initiate data simulation experiments, and find that ttreat results in fewer false discoveries compared to treat when using significance cut-offs of p =  or p = . but when the percentage of de genes in the simulated data is high, the advantage of ttreat over treat in terms of false discoveries is countered by ttreat missing more de genes. when the overall performance of treat and ttreat is compared by aucs, there is no difference on our simulated data. the treat bayesian approach shrinks the individual gene sample variances towards a pooled estimate, which works particularly well when the number of replicates is small; often no more than two replicates are available for microarrays. the nanostring counts in our dataset are not noisy and the throughput is much lower compared to a microarray experiment. so some of the gene-wise variances may be shrunk by too large a factor when applying treat on our nanostring data, as the hierarchical bayesian model is expecting a few outlier variances. second, when noisy data are expected, with many non-de genes showing signs of being increased or decreased, it is beneficial to use a safety margin around the chosen fc threshold in an analysis relative to a fc threshold. a good choice seems to set θ to τ +  or τ +  third, the particular choice of a fc threshold can influence the interpretation of biological data. the running fc model can be applied in order to avoid the subjective choice of a fc threshold, or when data vary strongly with expression levels. we used the same dataset to determine the fc thresholds for various gene expression levels and to apply ttreat, which may seem suboptimal. it would be ideal to determine fc thresholds on an independent dataset before using these thresholds in the subsequent analysis of the actual data. but often obtaining these independent data may be too expensive or impossible. we believe that performing the running fc model on the same dataset gives accurate results despite its data-driven nature in such circumstances. fourth, the choice of analytical tools must be driven by the research goals, the gene expression technique, and the hypothesis that is being tested. treat, ttreat, or ttreat either give a high number of false discoveries and few missed genes, only a few false discoveries for a higher number of missed genes, or an average but similar number of false discoveries and missed genes. the inevitable trade-off between false discoveries and missed genes must be evaluated dependent on the experimental context. in our case - nanostring studies of odorant receptor gene expression - we want to minimize the false discoveries, because the validation and follow-up experiments are laborious and expensive. but for other projects, the inverse may be desirable: the number of missed genes needs to be minimized, for instance in large-scale screening of molecules against a drug target, where missing a potential blockbuster drug cannot be afforded. in the running fc model, the choice of the fc percentile in each bin in the development of the lfc model as well as the choice of the k bins when applying the model affect the final results of the running fc model, thus the de genes discovered. depending on the percentage of de genes that is expected, percentile values ranging from p to p are deemed acceptable. regarding the number of k bins, a good coverage of the full expression range of the genes under analysis should be ensured. performance results on a single simulated dataset would have been subject to chance. repeating the simulations times permits an understanding of the variability in test performance that may be obtained between different simulations. the variability of our mean results is assessed by the 95% prediction intervals around the mean number of false discoveries on these data simulations. neither of the statistical tests shows outliers in performance that could contradict our conclusions. our data simulation procedure requires a “real” data input at the beginning only to get an idea about the distribution of gene error variances for nanostring. this input came from our data as described in the results. we also used an independent, published nanostring data set with ~ mouse genes to initiate the set of data simulation experiments again. we find that simulations initiated by these data lead to very similar results. we present approaches for testing the significance of gene expression relative to a fc threshold, with a focus on our nanostring data on or genes. when carrying out these hundreds of significance tests simultaneously, the issue of multiple testing can be raised. a significance analysis for microarrays test has been developed that encompasses both a test statistic for ranking and a way to control the false discovery rate and thus to adjust for multiple testing. the advantage of our approaches is that they provide p values that can be adjusted by a multiple testing method of choice, ranging from bonferroni to benjamini and hochberg. there is a wide variety of adjustment procedures for multiple testing. the choice of the most appropriate procedure for a given analysis is not straightforward. guidelines and criteria to be attributed to the multiple testing approaches have been described. in general, the formulation of a composite null-hypothesis forces treat, ttreat, and ttreat to rely on a more conservative way of p value calculation. the distribution of the resulting p values is thus skewed towards the larger values, which has consequences: if a method for adjusting multiple testing problems relies on the uniformity of the distribution of the p values, it cannot be used as this criterion is not met and we feel that when applied to nanostring data, these already conservative p values should not be adjusted by yet another conservative method. for the development of treat and initially also for t*, it was decided to use gene-wise linear models with arbitrary coefficients and contrasts of interest as a contextual and practical background; see also the limma package. the ttreat, ttreat and the running fc model are also based on gene-wise linear models. others have used linear models, or, more precisely, analysis of variance models to assess de genes in microarray experiments. the main difference is that all genes are included in a single linear model, presenting the advantage that the normalization process is combined with the data analysis and thus pre-normalization of the data is not necessary. obviously for such an anova model to be justified, its underlying model assumptions should be met, but they can be assessed quite straightforward. for a smaller nanostring codeset century, which consists of or genes and reference genes, we developed anova models that also dealt with assessing the de genes relative to a fc threshold. the negative controls could not be included in these anova models, as their residuals were making the general residual distribution heavily tailed and therefore non-gaussian. equality of variances was not met but within an acceptable borderline range on these genes. the results of these anova models were not satisfactory, as the genes that were identified as de genes were sometimes very different from the treat, ttreat and ttreat results, and less biologically meaningful. anova models including various genes in the same linear model can be applied on nanostring data and may be very useful for certain purposes. on our data, however, the gene-wise linear models and thus treat, ttreat, or ttreat were preferred over an approach that models all genes simultaneously.  our gene-wise statistical analysis of gene expression data with significance relative to a fc threshold gives reproducible and reliable results on nanostring data of odorant receptor genes, the largest gene family in the mouse genome. because it is difficult to set a biologically meaningful fc threshold in advance, we have developed methods that provide guidance in determining a stable fc threshold or in avoiding this choice altogether. by applying a two-stage approach, a safety margin around the fc threshold can be set, which is beneficial in certain situations. our running fc model identifies fc thresholds in an automated way, and is thus a more objective model leading to results that are more reproducible. this model is well suited for the problem of higher fc variances of lower expression values, thus avoiding a bias. 
statistical tests have been developed for microarray data to identify genes that are differentially expressed relative to a fold change threshold. gene-wise statistical analyses of gene expression data, for which the significance relative to a fold change threshold is important, give reproducible and reliable results on nanostring data of mouse odorant receptor genes. this set of tools may be useful for the analysis of other types of gene expression data. a statistical method, treat, has been developed for microarray data to assess formally if fold changes are significantly higher than a predefined threshold. a challenge in gene expression studies is the reliable identification of differentially expressed genes. gene expression studiesstatistical testdifferential gene expressionfold change criterionnanostring here we report that another approach, which we refer to as ttreat, is more appropriate for our nanostring data, where false discoveries lead to costly and time-consuming follow-up experiments. because it can be difficult to set in advance a fold change threshold that is meaningful for the available data, we developed methods that enable a better choice or avoid this choice altogether. we have recently applied the nanostring digital platform to study expression of mouse odorant receptor genes, which form with members the largest gene family in the mouse genome. methods that we refer to as ttreat and the running fold change model improve the performance of the statistical tests by protecting or selecting the fold change threshold more objectively. our objectives are, on these data, to decrease false discoveries when formally assessing the genes relative to a fold change threshold, and to provide a guided selection in the choice of this threshold.
 a challenge in gene expression studies is the reliable identification of differentially expressed genes. in many high-throughput studies, genes are accepted as differentially expressed only if they satisfy simultaneously a p value criterion and a fold change criterion. a statistical method, treat, has been developed for microarray data to assess formally if fold changes are significantly higher than a predefined threshold. we have recently applied the nanostring digital platform to study expression of mouse odorant receptor genes, which form with members the largest gene family in the mouse genome. our objectives are, on these data, to decrease false discoveries when formally assessing the genes relative to a fold change threshold, and to provide a guided selection in the choice of this threshold.  statistical tests have been developed for microarray data to identify genes that are differentially expressed relative to a fold change threshold. here we report that another approach, which we refer to as ttreat, is more appropriate for our nanostring data, where false discoveries lead to costly and time-consuming follow-up experiments. methods that we refer to as ttreat and the running fold change model improve the performance of the statistical tests by protecting or selecting the fold change threshold more objectively. we show the benefits on simulated and real data.  gene-wise statistical analyses of gene expression data, for which the significance relative to a fold change threshold is important, give reproducible and reliable results on nanostring data of mouse odorant receptor genes. because it can be difficult to set in advance a fold change threshold that is meaningful for the available data, we developed methods that enable a better choice or avoid this choice altogether. this set of tools may be useful for the analysis of other types of gene expression data. gene expression studiesstatistical testdifferential gene expressionfold change criterionnanostring 