Robust gene selection methods using weighting schemes for microarray data analysis

Background A common task in microarray data analysis is to identify informative genes that are differentially expressed between two different states. Owing to the high-dimensional nature of microarray data, identification of significant genes has been essential in analyzing the data. However, the performances of many gene selection techniques are highly dependent on the experimental conditions, such as the presence of measurement error or a limited number of sample replicates. Results We have proposed new filter-based gene selection techniques, by applying a simple modification to significance analysis of microarrays (SAM). To prove the effectiveness of the proposed method, we considered a series of synthetic datasets with different noise levels and sample sizes along with two real datasets. The following findings were made. First, our proposed methods outperform conventional methods for all simulation set-ups. In particular, our methods are much better when the given data are noisy and sample size is small. They showed relatively robust performance regardless of noise level and sample size, whereas the performance of SAM became significantly worse as the noise level became high or sample size decreased. When sufficient sample replicates were available, SAM and our methods showed similar performance. Finally, our proposed methods are competitive with traditional methods in classification tasks for microarrays. Conclusions The results of simulation study and real data analysis have demonstrated that our proposed methods are effective for detecting significant genes and classification tasks, especially when the given data are noisy or have few sample replicates. By employing weighting schemes, we can obtain robust and reliable results for microarray data analysis. Electronic supplementary material The online version of this article (10.1186/s12859-017-1810-x) contains supplementary material, which is available to authorized users.


Background
Microarray technologies allow us to measure the expression levels of thousands of genes simultaneously. Analysis on such high-throughput data is not new, but it is still useful for statistical testing, which is a crucial part of transcriptomic research. A common task in microarray data analysis is to detect genes that are differentially expressed between experimental conditions or biological phenotype. For example, this can involve a comparison of gene expression between treated and untreated samples, or normal and cancer tissue samples. Despite the rapid change of technology and the affordable cost for conducting whole-genome expression experiments, many past and recent studies still have relatively few sample replicates in each group, which makes it difficult to use typical statistical testing methods. These two problems, high dimensionality and small sample size problems, have triggered developments of feature selection in transcriptome data analysis [1][2][3][4][5][6][7][8][9]. These feature selection methods can be mainly classified into four categories depending on how they are combined with learning algorithms in classification tasks: filter, wrapper, embedded, and hybrid methods. For details and the corresponding examples of these methods, we refer the reader to several review papers [10][11][12][13][14][15][16][17][18]. As many researchers commented, filter methods have been dominant over the past decades due to its strong advantages, although they are the earliest in the literature [11-13, 15, 16]. They are preferred by biology and molecular domain experts as the results generated by feature ranking techniques are intuitive and easy to understand. Moreover, they are very efficient because they require short computation time. As they are independent of learning algorithms, they can give general solutions for any classifier [15]. They also have a better generalization property as the bias in the feature selection and that of the classifier are uncorrelated [19]. Inspired by its advantages, we focus on the filter method in this study.
One of the most widely used filter-based test methods is significance analysis of microarrays (SAM) [1]. It identifies genes with a statistically significant difference in expression between different groups by implementing gene-specific modified t-tests. In microarray experiments, some genes have small variance so their test statistics become large, even though the difference between the expression levels of two groups is small. SAM prevents those genes from being identified as statistically significant by adding a small positive constant to the denominator of the test statistic. This is a simple but powerful modification for detecting differentially expressed genes, considering the characteristics of microarray data. Since its establishment, the SAM program has been repeatedly updated. The latest version is 5.0 [20].
We also aim to develop methods for detecting significant genes based on a deeper understanding of microarray data. Even when researchers monitor an experimental process and control other factors that might have an influence on the experiment, biological or technical error can still arise in high-throughput experiments. For example, when one sample among a number of replicated samples gives an outlying result owing to a technical problem, variance of the gene expression becomes larger than expected and its test statistic becomes small. This is a major issue because it can lead to biologically informative genes failing to be identified as having a significant effect. Therefore, we here attempt to reduce this increase in variance for such cases by modifying the variance structure of SAM statistics, using two weighting schemes. It is also important to adjust the significance level of tests. Since we generally need to test thousands of genes simultaneously, the multiple testing problem arises. To resolve this problem, several methods have been suggested as replacements for the simple p-value; for example, we can use the family-wise error rate (FWER), false discovery rate (FDR) [1,21], and positive false discovery rate (pFDR) [22]. Among them, FDR, which is the expected proportion of false positives among all significant tests, is a popular method to adjust the significance level. It can be computed by permutation of the original dataset. The test procedures we propose in this paper also use FDR, the same as SAM.
Once a list of significant genes is established by a gene selection method, researchers may carry out further experiments such as real-time polymerase chain reaction to determine whether these reference genes are biologically meaningful. However, many genes may not be tested owing to limitations of time and resources. For example, even if hundreds of genes are included in a list of reference genes for a user-defined significance cutoff, researchers may just select a few top-ranked genes among them for further analyses. Therefore, it is very important that the genes are properly ranked in terms of their significance, especially for top-ranked genes [23,24]. As such, in this paper, we focus on improving test statistics for each gene and assessing how well each test method identifies significant genes.
For microarray data analysis, a comparison of the performance of gene selection methods is difficult because we generally do not know the "gold standard" reference genes in actual experiments. In other words, we do not know which genes are truly significant. This is a common problem encountered in transcriptome data analysis, so most studies have focused on comparing classification performances, which are determined by the combination of the feature selection and learning algorithm. As these results are clearly dependent on the performance of learning method, we cannot compare the effectiveness of feature selection techniques definitively [16]. Therefore, in this paper, we generate spike-in synthetic data that allow us to determine which genes are truly differentially expressed between two groups. For this, we suggest a data generation method based on the procedure proposed by Dembélé [25]. By performing such simulations, we can see how the performance changes depending on the characteristics of the dataset, such as sample size, the proportion of differentially expressed genes, and noise level. In this study, we focus on comparing performance according to noise level as our goal is to efficiently detect significant genes in a noisy dataset. To verify that our proposed methods can also compete with previous methods for actual microarray data, we use two sets of actual data that have a list of gold standard genes based on previous findings. All of these real datasets are publicly available and can be downloaded from a website [26] and R package [27]. In order to compare different gene selection methods, we also define two performance metrics that can be used when true differentially expressed genes are known. This paper is organized as follows. In the next section, we review the algorithm of SAM and propose statistical tests for microarray data that are modified versions of SAM, named MSAM1 and MSAM2. In addition, we explain our synthetic data generation method and suggest two performance metrics. In the results section, we describe our simulation studies and real data analysis. We compare SAM, MSAM1, and MSAM2 using 14 types of simulated dataset, which have different noise levels and sample sizes, and two sets of real microarray data. We next discuss the difference between the three methods in detail, focusing on FDR estimation. Additionally, we give the results of classification analysis using some topranked genes selected by each method. In the last section, we summarize and conclude this paper.

Methods
In this section, we briefly review the SAM algorithm [1] and propose new modified versions of SAM, focusing on calculating the test statistic.
SAM Let x ij and y ij be the expression levels of gene i in the jth replicate sample in states 1 and 2, respectively. For such a two-class case, the states of samples indicate different experimental conditions, such as control and treatment groups. Let n 1 and n 2 be the numbers of samples in these two groups, respectively. The SAM statistic proposed in [1] is defined as follows: where x i and y i are the mean expression of the ith gene for each group, x i ¼ P n 1 j¼1 x ij =n 1 and y i ¼ P n 2 j¼1 y ij =n 2 . The gene-specific scatter s i is defined as: where a = (1/n 1 + 1/n 2 )/(n 1 + n 2 − 2) and s 0 is a small positive constant called the fudge factor, which is chosen to minimize the coefficient of variation of d i . The computation of s 0 is explained in detail in [3]. Now let us consider the overall algorithm. The SAM algorithm proposed in [1] can be stated as follows: induced genes and d (down) as the largest d (i) among significantly suppressed genes. 9. The false discovery rate (FDR) is defined as the proportion of falsely significant genes among genes considered to be significant and can be estimated as follows: The algorithm consists of two parts: computation of the test statistic and determination of the cutoff for a given Δ. We will focus on the first of these parts and apply a simple modification to the computation of genespecific scatter s i to find a more robust test statistic. The numerator of the modified statistic and that of the original SAM statistic are the same. All of the procedures can be implemented using the samr package for Bioconductor in R. [20] described how to use the package and provided technical details of the SAM procedure.

Modified SAM
From one experiment [28], we observed several cases in which most of the results of gene expression are very close to each other, apart from one substantial outlier. As a result, the ranks of these genes from SAM are lower than expected. This prompted us to propose a new test method that has a different variance structure, leading to robustness on identifying informative genes in the presence of outliers. Throughout the paper, we use the term "outliers" to indicate "unusual observations".
Let us consider two cases with the following data: case 1: (5,5,5,5,8.54) and case 2: (3,4,5,6,7). For these two cases, the variance is the same, inferring that they have the same spread. However, even though the levels of variance are equal, in fact, we cannot say that the data points are similarly distributed. We believe that case 1 is more reliable than case 2. Our goal, therefore, is to propose a test statistic that has a more significant result for case 1 than for case 2. To minimize the effects of outliers among samples, we use the median instead of the mean and employ a weight function w when computing the test statistic, resulting in a less weight on an outlier sample that is far from other samples. A modified s i , s i , is defined as follows: Accordingly, our test statistic d~i is defined as follows: Methods modified by this approach might be particularly useful when detecting differentially expressed genes from noisy microarray data. The key idea is to reduce the impact of outliers when calculating the test statistic. We propose two different weight functions in this paper. The values of s i and d~i would differ quite markedly depending on the used weight function.

Modified SAM1 (Gaussian weighted SAM)
The weight function used in Modified SAM1 (MSAM1) is based on the Gaussian kernel, which is a widely used weight that decreases smoothly to 0 with increasing distance from the center. It is defined as follows: where ϕ is the probability density function of a stand- . The mean μ i is a gene-specific parameter such that μ i = median j (x ij ) and standard deviation σ is a data-dependent constant determined by the following procedure: first, m is defined as follows. m = max(|x ij − median j (x ij )|, |y ij − median j (y ij )|). It is calculated from given data. Second, p is a user-defined value between 0 and 1. Finally, given m and p, we can find the value of σ that satisfies the following equation: where F is the cumulative distribution function of a normal distribution. Therefore, m would approximately be the 100(1 − p)th percentile point of a normal distribution with mean 0 and standard deviation σ. As can be seen from Fig. 1, smaller p yields smaller σ. Therefore, smaller p makes the weight applied to outlier samples smaller.
On the other hand, as p increases, the results of original and modified SAMs become similar because the weight on the outlier is very similar to the weight on the nonoutliers. In this research, we set p = 0.001 since we found that this value is sufficiently small to reduce the effect of outliers.
For a better understanding of MSAM1, we here illustrate the weight function of MSAM1 and its application in detail. Let us consider Leukemia data [29]; for details of this data, see real data analysis section. The data consist of 38 samples (27 from ALL patients and 11 from AML patients) and 7129 genes. For simplicity and clarity, we randomly selected five samples for each sample type and applied SAM, MSAM1 with p = 0.01 and MSAM1 with p = 0.001. In order to compare weights given by each method, let us take one gene, M96326_rna1_at (Azurocidin). This gene would be a good example to clarify the difference between SAM and MSAM1 because it has an outlier sample. From Fig. 2, we can see that gene expressions in group 1 are  Table 1 and Fig. 3 show its gene expressions and weights computed by SAM and MSAM1. In Fig. 3, the lengths of 5 red dashed lines indicate the weights on the 5 observations. As we stated above, we can also see that smaller p makes the difference between weights applied to outlier and non-outlier samples greater.

Modified SAM2 (inverse distance weighted SAM)
This method uses Euclidean distance among the observations. The weight function used in Modified SAM2 (MSAM2) is defined as follows: is the Euclidean distance between the jth and kth samples of gene i. The reason that we use this weight function can be explained by the following example. Let us assume that there are 10,000 genes (i = 1 , 2 , … , 10000). Also, suppose there are 4 sample replicates (observations) in a group of the first gene (i = 1) and their gene expressions are x 11 , x 12 , x 13 and x 14 . Let w j be the weight on jth observation for j=1, 2, 3 and 4. In this case, the weights on these observations are as follows.
If x 11 , x 12 and x 13 are close to each other and x 14 is far from these 3 values, w 4 is much smaller than w 1 , w 2 and w 3 . Therefore, by using this weight function, we can give a smaller weight to an outlier. The further away an observation is from the others, the smaller weight is given.

Synthetic data generation
To run experiments, we need to generate synthetic gene expression data. These datasets should have characteristics similar to those of real microarray data to ensure that the results are reliable and valid. Two important characteristics of gene expression data, which are reported elsewhere [25,30,31] and also considered in this study, are as follows: 1. Under similar biological conditions, the level of gene expression varies around an average value. In rare cases, technical problems would result in values far away from this average. 2. Genes at low levels of expression have a low signalto-noise ratio.
The 'technical problems' mentioned in the first of these points are one possible explanation for outliers observed in microarray data. Since our goal is to develop methods that detect differentially expressed genes well in a noisy dataset containing outliers, we consider not only a dataset with little noise, but also a noisy dataset with outliers. We ensure that outliers are present at higher probability in several of the datasets to provide a wider range of comparisons among the different test methods. Basically, we follow the microarray data generation model by Dembélé [25], which uses a beta  distribution. In this article, we employ a beta and a normal distribution to generate data points, assuming that the levels of gene expression essentially follow such distributions. To allow outliers in generated data, we add a technical error term in our model; this term is mentioned in [25], but not used in their model. According to the noise level and distribution type, we consider four different simulation set-ups as follows: Scenario 1, noncontaminated beta; 2, contaminated beta; 3, noncontaminated normal; 4, contaminated normal. Therefore, data used in scenarios 1 and 3 have low noise level, and data used in scenarios 2 and 4 have high noise level. The step-by-step procedure for our data generation method is summarized as follows.
Step 1. Let n be the number of genes and n 1 and n 2 be control and treatment sample sizes, respectively.
Step 2. Generate z i from a beta (normal) distribution for i = 1 , 2 , … , n and transform the values, z i ¼ lb þ ub Âz i .
Step 3. For each z i , generate (n 1 + n 2 ) values as follows: Step 4. The final model is given by where the term s ij allows us to define differentially expressed genes. Their values are zero for the control group, s ij e N μ de ; σ 2 de À Á for genes with induced expression, and s ij e N −μ de ; σ 2 de À Á for genes with suppressed expression, where μ de ¼ μ min de þ Exp λ 2 ð Þ. n ij is an additive noise term, n ij e N 0; σ 2 n À Á . The final term t ij is used to define outlying samples by allowing non-zero values for some genes. The undefined parameters for each step can be set by the users. The values we use in this paper are as follows: λ 1 = 0.13, λ 2 = 2, μ min de ¼ 0:5 , σ de = 0.5, σ n = 0.4.
For these parameters, the influence of different parameter settings on the generated data is well explained elsewhere [25].

Scenario 1: Beta with low noise level
In this case, we generate data points from Beta(shape 1 , shape 2 ). shape 1 and shape 2 are two shape parameters of the beta distribution and we here set shape 1 = 2 and shape 2 = 4. We also set lb = 4, ub = 14. The values of t ij are zero for this case.

Scenario 2: Beta with high noise level
Here, we generate a noisier data than above data. The generation procedure is basically the same as the above case, except for allowing some non-zero t ij . To make outlying samples, we contaminate the data by adding gaussian noise to some treatment samples: For genes with induced or suppressed expression, t ij e N 0; σ 2 deo À Á for j = (n 1 + n 2 − n deo + 1) , … , (n 1 + n 2. ) where σ deo is a non-zero constant and n deo is the number of outlying samples. We here set σ deo = 1 and n deo = [0.2 × n 2 ] where [x] = m if m ≤ x < m + 1 for all integer m. For example, if there are five sample replicates in a treatment group, there can be one possible candidate as an outlier. Therefore, σ deo and n deo control the distribution and noise level of outlying samples. We believe that this set-up is reasonable because it does not destroy the original data structure while controlling the noise level of the data.

Scenario 3: Normal with low noise level
This scenario assumes that the levels of gene expression essentially follow a normal distribution, instead of a beta distribution. In this research, we use the normal distribution with mean 10 and standard deviation 1.5 for generated data points to be distributed between realistic bounds; the gene expression levels on a log2 scale after robust multichip analysis normalization usually vary between 0 and 20. We set lb = 0, ub = 1 in Step 2, which means that no transformation is applied.

Scenario 4: Normal with high noise level
To generate a noisier normal data, we use the same data generation procedure of Scenario 3, except for allowing some non-zero t ij in Step 4. The structure of t ij is the same as in Scenario 2.

Performance metrics
To compare the performance of several methods, we need several evaluation measures. Since we know which genes are differentially expressed in our simulated datasets, we can define two performance metrics as follows, measuring how well each method identifies these TRUE genes. Prior to define metrics, let G up ={i: gene i the expression of which is truly significantly induced} and G down ={i: gene i the expression of which is truly significantly suppressed}.

Rank sum (RS)
We define the rank sum (RS) of TRUE genes as follows: where I(•) is an indicator function. The reason for determining the ranks of genes with high and low expression is that the SAM procedure uses such a method when detecting genes of the two groups. We use the absolute value of test statistics because test statistics of genes with suppressed expression have negative values. For RS, lower values indicate better performance.

Top-ranked frequency (TRF)
The top-ranked frequency (TRF) of TRUE genes is computed by Here, r denotes the rank cutoff and is set to be smaller than the number of observations in G up and G down . For a given cutoff r, TRF computes the number of TRUE genes ranked within r. For TRF, higher values indicate better performance.
To understand the performance metrics better, let us consider the following case. We have 100 genes and 10 TRUE genes among them. Assume that we obtain a topranked gene list as shown in Table 2

Simulation studies
In this section, we compare gene selection methods using synthetic datasets. We consider four scenarios described above. For each scenario, we consider 7 different combinations of n 1 and n 2 in order to take into account the affects of sample size and class imbalance on gene selection performance as follows: (n 1 , n 2 ) = (5, 5) , (5, 10) , (10, 5) , (10, 10) , (10,15) , (15,10) and (15,15). For all scenarios, we assume that there are 2% target genes (1% up-regulated and 1% down-regulated genes) among the total of 10,000 genes. For simplicity, let us assume that the first 100 genes are downregulated and last 100 genes are upregulated. Then, we can describe the structure of our simulation data as shown in Fig. 4. This example illustrates the structure of noisy data containing outliers. In this case, the last two samples are outlying samples among 10 treatment samples of 200 target genes. There are five different distributions of data points: A, B, C, D, and E. For 9800 nontarget genes, the distributions of the control and treatment samples are the same (A). The first 100 downregulated genes are generated from two distributions (B and C) and the last 100 upregulated genes are also generated from two distributions (D and E). Groups C and E indicate outlier samples. If there are no outliers in the dataset, B is equivalent to C and D is equivalent to E. The empirical density plot of each group is shown in Fig. 5. For visualization, we use 5000 data points to ensure equivalent density of the points for each group (A, B, and C), that is, with a 1:1:1 ratio, not using the original ratio among the three groups. We conduct simulation studies using synthetic data and compare the results using three metrics; two of them are RS and TRF, which were defined above, and the third is AUC. AUC is the area under a receiver operating characteristic (ROC) curve. Therefore, this value falls between 0 and 1, and higher values indicate better performance. We consider five gene selection methods, named SAM, SAM-wilcoxon, SAM-tbor, MSAM1 and MSAM2. SAM-wilcoxon is the Wilcoxon version of SAM [20,32]. SAM-tbor is basically the same with SAM, except for applying a simple trim-based outlier removing algorithm to data prior to running SAM. In this study, we remove the largest and smallest observations from each sample type. Figs. 6 and 7 display the average Fig. 4 An example of simulated data structure. Each row and each column of this data frame correspond to a gene and a replicate sample, respectively, so we have a 10,000 × 20 data matrix in this study. We assume that there are 2% target genes (1% up-regulated and 1% down-regulated genes) among the total of 10,000 genes, and ten replicates in each group. There are five different distributions of data points: A, B, C, D, and E; groups C and E indicate outlier samples  performance of 100 simulations for each method on the three metrics. Table 3 shows numerical results of 4 cases. The best performance on each metric is shown in boldface. In scenario 1, the original SAM always outperform SAM-wilcoxon and SAM-tbor. Although SAMtbor show better performance than SAM in some cases of scenario 2, its performance is worse than those of MSAMs. As can be seen from the figures and table, our proposed methods show better performance than three versions of SAM in all cases. In particular, modified SAMs are much better when given data is noisy (scenario 2, compared to scenario 1) and is a little better for less noisy cases. We can also see that our methods show more robust performance in all cases. When there is two outliers among ten samples, the number of target genes found by original SAM is reduced by 2-17%, whereas that found by MSAMs is reduced by 1-8%. In particular, when n 1 = 5, n 1 = 10 in scenario 2, SAM fails to detect 90 genes among the 200 TRUE genes, whereas MSAM2 fails to detect only 60 genes on average. Simulation results of scenarios 3 and 4 are in Additional file 1. These results are very similar with those of scenarios 1 and 2; MSAMs always perform better than three versions of SAM.
In real data analysis sections, we only consider SAM, MSAM1, and MSAM2, all of which show good performance in simulation studies; we found that SAM-wilcoxon and SAM-tbor are worse than the original SAM in the previous section. Moreover, we cannot apply SAM-tbor to this data because this data has only three sample replicates in each group. Like this case, we can see that such a trim-based method is limited in its applications. Tables 4 and 5 show the rank of 11 reference genes that are differentially expressed between the control group and the treatment groups (dtri6 and dtri10, respectively). The last row in each table indicates the rank sum of these 11 genes. As we can see, MSAM2 shows the best performance because the rank sum of this method is the smallest among those of the three gene Note: the best performance in each case is shown in bold type selection methods. In particular, MSAMs improve the rank of the genes named fgd4-170_at and fgd159-500_at. For each of these genes, the result for one of their treatment samples is far from those for the other two samples. From the analysis, it can be asserted that our proposed methods efficiently identify the genes whose replicate samples contain an outlier, such as fgd4-170_at and fgd159-500_at.

Real data analysis 2: Leukemia
Leukemia is a cancer of the bone marrow, where blood cells are made. In leukemia, abnormal blood cells are produced in the bone marrow and crowd out other normal blood cells. Depending on the type of abnormal blood cells that are multiplying, leukemia can be classified as acute lymphocytic leukemia (ALL) or acute myeloid leukemia (AML). Identifying the type of leukemia is very important because patients should receive different treatments according to the disease type.
[29] studied a generic approach to cancer classification based on gene expression and provided a list of 50 significant genes for classifying ALL and AML. After this study, this dataset has been widely used in transcriptomic analysis, e.g., [34,35]. This data are available in the golubEsets library in Bioconductor [27]. The original data consist of 38 samples (27 from ALL patients and 11 from AML patients) and 7129 genes. We randomly selected five, seven, and ten samples for each sample type and repeated this experiment 100 times for averaging because biological experiments usually have a small number of samples owing to limitations of time and resources. It is thus important that a method shows good performance even if the sample size is small. The simulation results are shown in Table 6. In this table, RS and TRF values of three gene selection methods, which were computed by using 50 genes that are considered informative in [29] over 100 trials. For each case, the best performance is shown in boldface in the table. As we can see, MSAM1 or MSAM2 performs better than SAM in terms of RS and TRF, regardless of rank cutoff values. The overall performance of SAM and MSAM1 are very similar, but MSAM1 always performs slightly better than SAM. In the point of view of sample size, MSAM2 outperform SAM and MSAM1 when the sample size is very small, e.g., 5, and MSAM1 performs better than SAM and MSAM2 when the sample size is moderate, e.g., 7 and 10. As the sample size increases, all of the three methods identify informative genes better.

FDR comparison
In this section, we discuss the FDR estimation procedures of SAM, MSAM1, and MSAM2. FDR is used in SAM procedure in order to deal with a multiple testing problem. The SAM interface in R, samr package [20], provides a significant gene list based on the FDR value that is estimated by its internal function. We also construct our own interface for MSAMs in R, based on the samr package, in order to allow for users to apply our proposed methods to their transcriptome research; see Additional file 2. Users start the procedure by setting their desired FDR value (for example, 0.2). We will call this value 'estimated FDR'. Based on the estimated FDR, our procedure calculates the value of corresponding Δ and identifies potentially significant genes. In real applications, we do not know TRUE FDR, so the estimated FDR is used as a substitute for TRUE FDR. If the estimated value is different from the true value, the number of genes that are detected using the estimated FDR is larger or smaller than the true number. Therefore, users may be interested in how well SAM and MSAMs  Since we know the number of TRUE significant genes in our simulated datasets, we can compare the estimated FDR and TRUE FDR in simulation study. After 100 simulations, we draw a scatter plot of the TRUE FDR versus the estimated FDR by calculating the average values of the TRUE FDR for each estimated FDR. We next draw a smooth curve close to the scatter plot for scenarios 1 and 2 to find the estimation accuracy at various levels of FDR. In particular, the estimation accuracy at low FDR is important since researchers generally set FDR at a small value so as to avoid having a large proportion of falsely significant genes among the detected genes. For this reason, we only show the results when the estimated FDR is lower than 0.5. Figure 8 displays the results; see the top two plots. As we can see, SAM estimates the TRUE FDR very accurately and two modified SAMs slightly overestimate the TRUE FDR. In other words, our methods have conservative property in their FDR estimation. However, the conservative estimation of FDR may not cause serious problems for the analysis when we use FDR as an upper bound of a tolerable error [36].
For such an analysis, the more important thing is how many non-significant genes are included in the detected genes. Because the truths are known in the simulated data, we can calculate the number of falsely detected genes among the identified genes. With the same number of total positives, the method with the smallest number of false positives is the best [36]. Using the plotting method described above, a smooth curve of the number of false positive genes versus the total number of identified genes are drawn. Figure 8 shows the results From the figure, we can see that MSAM1 and MSAM2 gives smaller number of false positive genes than SAM across all noise level and the total number of identified genes. From the results, we can say that MSAMs are better than SAM because they includes the less number of false genes in the selected gene subset.
When we estimate FDR, we calculate both median FDR and mean FDR to determine which estimate more closely approximates the true value. Since the original samr interface provides the median FDR and 90th percentile FDR only, we modified its estimation function and obtained the median and mean values of FDR. As a result, we found that the median FDR was closer than the mean FDR to the TRUE FDR for all methods. This coincides with results published elsewhere [37], in which the median FDR was recommended as a criterion for gene selection methods when the estimated proportion of differentially expressed genes is greater than 1%, regardless of the sample size. Based on these results, we use the median value instead of the mean value when estimating FDR.

Classification analysis
Once important genes are identified from thousands of genes, they can be used to predict two different experimental states or responses (for example, cancer and normal). Therefore, we also examine how well a few top genes selected by each method identify the true classes. We attach these results in Additional file 3. In this file, we introduce 4 datasets we used and explain the construction of classifiers, 6 gene selection methods, 3 performance metrics to be considered in this study. Our comments on the results are also included. As can be seen in the file, our proposed methods, MSAMs, show quite good performances in all cases. In this additional Note: the best performance for each rank cutoff is shown in bold type section, we prove their competitiveness in classification tasks, not only in gene selection tasks.

Discussion
In transcriptome data analysis, most studies have been devoted to developing filter-based methods that are the simplest and fastest, and most computationally efficient. Hybrid methods, which are generally the combination of filter and wrapper methods, have recently gained popularity in the literature [13]. These methods consist of two steps: First, relevant features are selected by a filter method and the remaining features are eliminated. Second, a wrapper method verify these features and determine the final feature set that gives high classification accuracy [16]. In this point of view, filter methods have a lot of flexibility as they can be combined with not only any learning algorithm, but also any gene selection method, such as a wrapper method, resulting in a hybrid method. The performance of a hybrid method relies totally on the combination of filter and wrapper methods as well as the classifier [18]. We believe that accurate gene selection by filter methods clearly allow better classification accuracy. Therefore, our new filter-based methods will be useful not only in gene selection, but also constructing a good classifier in microarray applications. Our experiments showed the efficiency of our methods; it was demonstrated that when the same number of genes were selected, our methods included the less number of false genes than the conventional method. Our results also strongly suggest that these newly proposed methods outperform the conventional method and show quite consistent performance, even with a high noise level and a small sample size. Given that noisy data and a small sample size are commonly encountered in microarray studies [30,[38][39][40], we believe that our methods will prove useful.
This research was based on the existing interface of SAM that was modified to apply our proposed methods. This modified version of the samr package is available in Additional file 2. We attempted to find a balance between flexibility and control in the usage of our methods by allowing users to set particular parameters and by minimizing the number of modifications to the original