Abstract
Background
One important application of microarray experiments is to identify differentially expressed genes. Often, small and negative expression levels were clippedoff to be equal to an arbitrarily chosen cutoff value before a statistical test is carried out. Then, there are two types of data: truncated values and original observations. The truncated values are not just another point on the continuum of possible values and, therefore, it is appropriate to combine two statistical tests in a twopart model rather than using standard statistical methods. A similar situation occurs when DNA methylation data are investigated. In that case, there are null values (undetectable methylation) and observed positive values. For these data, we propose a twopart permutation test.
Results
The proposed permutation test leads to smaller pvalues in comparison to the original twopart test. We found this for both DNA methylation data and microarray data. With a simulation study we confirmed this result and could show that the twopart permutation test is, on average, more powerful. The new test also reduces, without any loss of power, to a standard test when there are no null or truncated values.
Conclusion
The twopart permutation test can be used in routine analyses since it reduces to a standard test when there are positive values only. Further advantages of the new test are that it opens the possibility to use other test statistics to construct the twopart test and that it avoids the use of any asymptotic distribution. The latter advantage is particularly important for the analysis of microarrays since sample sizes are usually small.
Background
The addition of a methyl group at the carbon5 position of cytosine is a modification of DNA called DNA methylation. In mammalian cells, DNA methylation is essential for proper development [1]. The methylation patterns of tumor cells are altered compared to those of normal cells, moreover, there are also differences between different types of cancer as shown for subtypes of leukemia [2] and lung cancer [3]. Thus, DNA methylation analysis promises to become a powerful tool in cancer diagnosis [4].
DNA methylation data can be obtained using the MethyLight technology [5]. When the tested region is not or only partially methylated the result is negative (undetectable methylation, null values). In contrast, samples that show methylation will have a value greater than 0 [4]. Thus, DNA methylation data obtained with MethyLight have a clump of zero observations and a continuous nonzero part. For such a data structure, twopart models as proposed by Lachenbruch [68] are applicable. In that approach, the test statistic is the sum of two squared statistics, one comparing the proportions of zeros and one comparing the positive values. For example, one can use the binomial test and the Wilcoxon rank sum test. The asymptotic null distribution of the sum of the squares of the two test statistics is χ^{2 }with two degrees of freedom (df = 2).
In microarray data it is relatively common that small and negative expression levels were clippedoff to be equal to an arbitrarily chosen cutoff value. Recent examples used different cutoff values: 1, 20, and 50, respectively [911]. Aside from the fact that negative values make no biological sense, there are, with regard to oligonucleotide arrays from Affymetrix, two primary reasons for truncating the values [12]. First, spots at the low intensity range are generally more vulnerable to noise, thus, it is thought that the technology produces a poor discrimination at low levels of expression [13]. Second, the focus is on expression of identified genes and expressed sequence tags. Differences at negative or low values may result from differences in binding to the mismatch probes. Since it is generally not known what binds to the mismatch probes, the differences at negative or low values cannot be attributed to target genes.
Due to the truncation there are two different types of data: truncated values and original observations. Since the truncated values are not just another point on the continuum of possible values, it would be inappropriate to use a standard statistical method that would treat all values equally [14]. The two different types of data should be analyzed separately. Therefore, the twopart model, comparing the proportions of truncated values and the distribution of positive values, is applicable. Note that negative expression levels are not possible when the Affymetrix Microarray Suite (MAS) 5.0 software is used. However, small values are possible and may be truncated.
Since, in microarray experiments, the sample sizes, i.e. the numbers of replications, are usually very small [15,16], the use of the asymptotic distribution of Lachenbruch's twopart test statistic may be questionable. Thus, we carried out permutation tests with the twopart statistic. For a permutation test all possible permutations under the null hypothesis are generated. In our situation we permute the group labels for the whole sample, i.e. for truncated values (or the null values in case of methylation data) and original observations. Then, the test statistic is calculated for each permutation. The null hypothesis can then be accepted or rejected using the permutation distribution of the test statistic, the pvalue being the probability of the permutations giving a value of the test statistic as supportive or more supportive of the alternative than the observed value [17,18]. Thus, inference is based upon how extreme the observed test statistic is relative to other values that could have been obtained under the null hypothesis.
We found that, in the case of a twopart model, the permutation test is not only a way to avoid the use of an asymptotic distribution, but also is a more powerful test, i.e. a test that produces, on average, smaller pvalues. In addition, the permutation test reduces, without any loss of power, to a single test if no truncated (or null) values were present. Thus, the proposed test is applicable in routine use whether or not truncated (or null) values occur. After the definition of the tests in the following section, we present our findings for DNA methylation data and microarray data. We then confirm the results using simulations.
Twopart tests
As briefly mentioned above a twopart test statistic is the sum of two squared statistics, one comparing the proportions of truncated values and one comparing the positive values. Let n_{1 }and n_{2 }be the numbers of independent observations regarding one gene (or one region in case of methylation data, respectively), for two groups to be compared. The observed numbers of truncated values (or null values in case of methylation data) in the two groups are denoted by m_{1 }and m_{2}. To compare these numbers m_{1 }and m_{2 }Lachenbruch [6] used the statistic
where , , and . Under the null hypothesis the proportions of truncated values are not different between the two groups, and B^{2 }is asymptotically χ^{2}distributed with df = 1. B^{2 }is always well defined, unless there are only truncated values in both groups or no truncated values at all. For these two extreme cases we set B^{2 }= 0.
For the second part in Lachenbruch's twopart model one can use different tests, Lachenbruch [68] considered the Wilcoxon rank sum test, Student's t test, and the KolmogorovSmirnov test. The use of the latter test in a twopart model, however, was too liberal, the type I error rate was close to 0.065 (for a significance level of α = 0.05 [7,8]). Both other tests work well. Here, we apply the Wilcoxon test for two reasons. This test was used in the original analyses of the data we use below [3,11]. Nonparametric tests based on ranks are more appropriate for nonnormally distributed data such as microarray data [19].
The standardized rank sum statistic based on the nontruncated (or positive values in case of methylation data) values is defined as
where RS is the rank sum, i.e. the sum of the ranks in group 1. When there are ties within the nontruncated values the denominator slightly changes [[20], p. 109]. For the extreme case that there are only truncated values in at least one group we set W = 0.
The test statistic for the twopart test is X^{2 }= B^{2 }+ W^{2}. Under the null hypothesis of no difference between the groups, X^{2 }is asymptotically χ^{2}distributed with df = 2 [6,7]. Alternatively, a permutation test can be performed with the statistic X^{2}. This test, called twopart permutation test here, is a permutation test based on the sum statistic X^{2}. It is carried out by permuting the group labels for the whole sample. Thus, all observations, truncated and nontruncated values (or null and positive values in case of methylation data) are reallocated to the groups. When performing this twopart permutation test the exact permutation distribution of X^{2 }is determined. This distribution, computed by generating all possible permutations or, for the pvalues given below in Table 1, using a simple random sample of 20,000 permutations, is used to compute the pvalue. Since it is a permutation test based on the sum X^{2}, it is neither necessary to determine the permutation distributions of the summands B^{2 }and W^{2 }nor to calculate the pvalues of the univariate tests related to B^{2 }and W^{2}.
Table 1. pvalues of the original twopart test and the twopart permutation test for seven regions in lung cancer cell lines; DNA methylation data from Siegmund et al. [4]
Application to actual methylation data
We use DNA methylation data from 7 regions and 87 lung cancer cell lines, 41 lines are from small cell lung cancer and 46 lines from nonsmall cell lung cancer [3,4]. The proportion of positive values for the different regions ranges from 39 to 100% for the small cell lung cancer and from 65 to 98% for the nonsmall cell lung cancer. The data are available at http://wwwrcf.usc.edu/~kims/SupplementaryInfo.html webcite. Siegmund et al. [4] transformed the data by standardizing the positive values on the naturallog scale. However, for the tests applied here, this transformation has no influence.
Table 1 presents the pvalues of the two tests. For most regions the twopart permutation test gives a smaller pvalue than the original twopart test. The only exception is the region APC. However, the original twopart test's pvalue for this region is 0.1684 and, for a pvalue of this size, a small change in the value is usually of no importance.
Application to actual microarray data
Tschentscher et al. [11] performed an experiment with HGU95Av2 oligonucleotide arrays in order to compare patients with uveal melanomas with and without monosomy 3. Expression values were calculated by use of the MAS 4.0 software. The data are available at http://www.uniessen.de/humangenetik/download webcite. The sample size in this microarray experiment is 10 per group. As mentioned above, expression levels below 50 were set to 50. This data truncation occurred in 2,215 (28%) out of 7,902 genes. First, we consider these 2,215 genes. Figure 1 displays the number of genes for the different number of truncated values per gene.
Figure 1. Number of genes with the given number of truncations per gene (data from the microarray experiment of Tschentscher et al. [11], only genes with at least one truncation)
Table 2 shows the frequencies of different size groups of the pvalues. Often, the twopart permutation test gives a smaller pvalue than the original twopart test. For instance, for 19 genes the pvalue is ≤ 0.001 when the latter test is applied. The permutation test gives a pvalue ≤ 0.001 for these 19 and for 39 additional genes. As usual, the majority of genes do not show any indication of being differentially expressed. For these genes a slight change in the pvalue is of no importance. Thus, in Table 3, we consider the genes for which the pvalue of the original twopart test is ≤ 0.1. Out of the 2,215 genes 514 remain. As shown in Table 3 the pvalues of the twopart permutation test are, on average, distinctly smaller than those of the original twopart test.
Table 2. Frequencies of different size groups of the pvalues of the original twopart test and the twopart permutation test for genes with at least one truncation; data from the microarray experiment of Tschentscher et al. [11]
Table 3. Differences between the pvalues of the original twopart test and the twopart permutation test for genes with at least one truncation and small pvalues (i.e. the pvalue of the original test must be as small as mentioned under condition), a positive difference means that the twopart permutation test has a smaller pvalue than the original test; data from the microarray experiment of Tschentscher et al. [11]
For a large proportion of genes (72% in this data set) there are no truncated values. In that case, the twopart statistic reduces to the sum 0 plus the squared standardized rank sum W^{2}. Of course, one has to define a priori whether the twopart test or the Wilcoxon test will be used to analyze the data. If the original twosided (asymptotic) Wilcoxon test were chosen and applied, one could compare W^{2 }with critical values from a χ^{2 }distribution with one degree of freedom (df). If the twopart were chosen there is df = 2 and, in case of no truncated values, power is lost compared to the original Wilcoxon test. For instance, the 95% percentile of the χ^{2 }distribution with df = 1 is 3.84, but it is 5.99 for df = 2. The permutation test using the sum statistic X^{2 }does not suffer from this power loss: When there are no truncated values there is, of course, no difference in the proportions of truncated values, and the test statistic X^{2 }is, for every permutation, the sum 0 + W^{2}. Thus, the twopart permutation test reduces to the exact twosided Wilcoxon rank sum test when there are no truncated values. Consequently, this permutation test does not only give smaller pvalues, but it is also applicable in routine use whether or not truncated values are present.
Simulation study
The two different tests, the original twopart test and the twopart permutation test, were compared in a Monte Carlo simulation study performed using SAS version 8.2, 5,000 simulation runs were generated for each configuration. The sample size of 10 per group was chosen as in the microarray experiment presented above. In some configurations some randomly chosen values were set to 0 according to binomial distributions with the probabilities p_{1 }and p_{2}, and the remaining observations were generated according to a lognormal distribution (with median 1 and σ = 1). Then, the values of one group were shifted if applicable. In some other configurations all observations were generated according to the lognormal distribution and, in one group, shifted. Then, values smaller than a cutoff value were truncated.
The type I error rates of the two tests are very similar. With e.g. p_{1 }= p_{2 }= 0.3 and a significance level of α = 0.05 the simulated type I error rates were 0.049 for both the original twopart test and the twopart permutation test.
Table 4 displays results for situations with a difference between the two groups. As above, only those comparisons were regarded for which there is some indication of a difference, i.e. a pvalue ≤ 0.1 of the original twopart test. In all considered configurations the median of the difference between the pvalues of the original twopart test and the twopart permutation test is positive. The finding that the pvalues of the permutation test are smaller corresponds to a higher power of this test. The power is given in Table 5, as shown the power of the twopart permutation test is at least as high as that of the original twopart test. There is only one exception, the latter test is slightly more powerful in one situation, i.e. when the proportion of zeros is higher in group 1 and the positive values are larger in group 1.
Table 4. Results of the simulation study: Differences between the pvalues of the original twopart test and the twopart permutation test for data sets with small pvalues (i.e. pvalue of the original test ≤ 0.1), a positive difference means that the twopart permutation test has a smaller pvalue than the original test; p_{1 }and p_{2 }are the probabilities for zero values and the positive values in group 1 are shifted by μ; 5,000 data sets with n_{1 }= n_{2 }= 10 were generated for each configuration
Table 5. Results of the simulation study: Powers of the original twopart test and the twopart permutation test; 5,000 data sets with n_{1 }= n_{2 }= 10 were generated for each configuration (notation as in Table 4, significance level α = 0.05)
Discussion
Previous research demonstrated that a twopart test is appropriate and powerful in the presence of a clump of zero observations (i.e. truncated or null values). In this paper we propose a permutation test for such situations with two types of data. Usually, nonparametric tests can be performed based on an asymptotic distribution or based on a permutation null distribution. The two approaches often give similar results, especially when the sample sizes are large. However, in the case of a twopart test one cannot simply replace the asymptotic distributions of B^{2 }and W^{2 }by the exact permutation distributions. If so, one would compute two exact pvalues although the aim of a twopart test is to receive one pvalue that combines information from both parts. Therefore, the twopart permutation test uses the exact permutation distribution of the sum statistic X^{2}. That this permutation distribution of the sum is generated rather than to simply replace the asymptotic distributions of the summands by their exact permutation distributions may be the reason why the permutation test is more powerful.
A disadvantage of a permutation test is that it can be computerintensive. However, this issue is less relevant now due to faster algorithms [[18], chap. 13] and the advent of highspeed PCs. Furthermore, one can carry out a permutation test based on a random sample out of the possible permutations, as we did for the DNA methylation data (see Table 1).
In microarray experiments it is common to investigate thousands of genes simultaneously. The approach presented here for the identification of differentially expressed genes is to consider a univariate testing problem for each gene. A correction for the multiplicity of genes is a subsequent step, that is, like the previous step of normalizing the data, outside the scope of this paper. A common approach to the multiplicity problem is to consider a procedure for testing the genes simultaneously for differential expression with the test on an individual gene being implied in the simultaneous test. For such a procedure different proposals have been made recently. For instance, there are methods based on the pvalues of the tests from individual genes [2123]. In a similar manner, the multiplicity of regions can be managed in DNA methylation data.
Conclusion
Aside from the shown improvement in power, the proposed twopart permutation test has three important advantages. First, it avoids the use of any asymptotic distribution and, therefore, can safely be applied in case of small sample sizes that are common in microarray experiments. Second, it reduces without any loss of power to the exact Wilcoxon test if there were no truncated (or zero) values. Thus, it can be used in routine analyses. Third, the permutation test opens the possibility to use other tests to construct the twopart test. Thus, tests with unknown or nonstandard null distributions can be used. For instance, one could replace the Wilcoxon test by the BaumgartnerWeißSchindler test [24] that was recently recommended for the analysis of gene expression data [19].
Authors' contributions
MN performed the statistical analyses and drafted the manuscript. TB prepared the microarray data. TB and KHJ participated in the design of the simulation study and helped to draft the manuscript. All authors read and approved the final manuscript.
References

Tsou JA, Hagen JA, Carpenter CL, LairdOffringa IA: DNA methylation analysis: a powerful new tool for lung cancer diagnosis.
Oncogene 2002, 21:54505461. PubMed Abstract  Publisher Full Text

Model F, Adorjan P, Olek A, Piepenbrock C: Feature selection for DNA methylation based cancer classification.
Bioinformatics 2001, 17(Suppl 1):S157S164. PubMed Abstract  Publisher Full Text

Virmani AK, Tsou JA, Siegmund KD, Shen LYC, Long TI, Laird PW, Gazdar AF, LairdOffringa IA: Hierarchical clustering of lungcancer cell lines using DNA methylation markers.
Cancer Epidemiology, Biomarkers & Prevention 2002, 11:291297. PubMed Abstract  Publisher Full Text

Siegmund KD, Laird PW, LairdOffringa IA: A comparison of cluster analysis methods using DNA methylation data.
Bioinformatics 2004, 20:18961904. PubMed Abstract  Publisher Full Text

Eads CA, Danenberg KD, Kawakami K, Saltz LB, Blake C, Shibata D, Danenberg PV, Laird PW: MethyLight: a highthroughput assay to measure DNA methylation.
Nucleic Acids Research 2000, 28:E32. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lachenbruch PA: Comparison of twopart models with competitors.
Statistics in Medicine 2001, 20:12151234. PubMed Abstract  Publisher Full Text

Lachenbruch PA: Analysis of data with excess zeros.
Statistical Methods in Medical Research 2002, 11:297302. PubMed Abstract  Publisher Full Text

Jelinek DF, Tschumper RC, Stolovitzky GA, Iturria SJ, Tu Y, Lepre J, Shah N, Kay NE: Identification of a global gene expression signature of Bchronic lymphocytic leukemia.
Molecular Cancer Research 2003, 1:346361. PubMed Abstract  Publisher Full Text

Küppers R, Klein U, Schwering I, Distler V, Bräuninger A, Cattoretti G, Tu Y, Stolovitzky GA, Califano A, Hansmann ML, DallaFavera R: Identification of Hodgkin and ReedSternberg cellspecific genes by gene expression profiling.
Journal of Clinical Investigation 2003, 111:529537. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tschentscher F, Hüsing J, Hölter T, Kruse E, Dresen IG, Jöckel KH, Anastassiou G, Schilling H, Bornfeld N, Horsthemke B, Lohmann DR, Zeschnigk M: Tumor classification based on gene expression profiling shows that uveal melanomas with and without monosomy 3 represent two distinct entities.
Cancer Research 2003, 63:25782584. PubMed Abstract  Publisher Full Text

Ibrahim JG, Chen MH, Gray RJ: Bayesian models for gene expression with DNA microarray data.
J Am Stat Assoc 2002, 97:8899. Publisher Full Text

Tamayo P, Slonim D, Mesirov J, Zhu Q, Kitareewan S, Dmitrovsky E: Interpreting patterns of gene expression with selforganizing maps: methods and application to hematopoietic differentiation.
Proc Nat Acad Sci USA 1999, 96:29072912. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Delucchi KL, Bostrom A: Methods for analysis of skewed data distributions in psychiatric clinical studies: working with many zero values.
American Journal of Psychiatry 2004, 161:11591168. PubMed Abstract  Publisher Full Text

Gadbury GL, Page GP, Heo M, Mountz JD, Allison DB: Randomization tests for small samples: an application for genetic expression data.

Zhao Y, Pan W: Modified nonparametric approaches to detecting differentially expressed genes in replicated microarray experiments.
Bioinformatics 2003, 19:10461054. PubMed Abstract  Publisher Full Text

Manly BFJ: Randomization, Bootstrap and Monte Carlo Methods in Biology. 2nd edition. Chapman and Hall, London, U.K; 1997.

Good PI: Permutation Tests. 2nd edition. Springer, New York, NY; 2000.

Neuhäuser M, Senske R: The BaumgartnerWeißSchindler test for the detection of differentially expressed genes in replicated microarray experiments.
Bioinformatics 2004, 20:35533564. PubMed Abstract  Publisher Full Text

Hollander M, Wolfe DA: Nonparametric statistical methods. 2nd edition. Wiley, New York, NY; 1999.

Zaykin DV, Zhivotovsky LA, Westfall PH, Weir BS: Truncated product method for combining Pvalues.
Genetic Epidemiology 2002, 22:170185. PubMed Abstract  Publisher Full Text

Dudbridge F, Koeleman BPC: Rank truncated product of Pvalues, with application to genomewide association scans.
Genetic Epidemiology 2003, 25:360366. PubMed Abstract  Publisher Full Text

Storey JD, Tibshirani R: Statistical significance for genomewide studies.
Proceedings of the National Academy of Sciences USA 2003, 100:94409445. Publisher Full Text

Baumgartner W, Weiß P, Schindler H: A nonparametric test for the general twosample problem.