Liu et al. [1] introduce an interesting alternative methodology for clustering genes using order-restricted inference methodology. Their strategy differs from that of [2,3] in that they avoid the bootstrap computation of p-value for identifying significant genes. Instead, they use an information theoretic model selection criterion to assign genes to temporal patterns (clusters). Once genes are clustered, they evaluate the reliability of each cluster using a bootstrap algorithm along the lines of Kerr and Churchill [5]. Liu et al. [1] make several statements about our methodology [2,3] that are erroneous and require clarification.

Corresponding author: Baoxue Zhang

We are delighted to see a thorough correspondence to our paper [7]from Dr. Peddada and his coauthors. In the correspondence, Peddada et al. made several comments about our paper on the comparison between their clustering approach, ORIOGEN [8]and our ORICC algorithm in [7].

The first comment is that the computation time reported in [7]about ORIOGEN is not accurate, and it could be we either misinterpreted details of ORIOGEN or our coding is extremely inefficient.

We first need to point out that, in our paper, we compare between ORICC and Peddada's method that refers to the algorithm in [9]. ORIOGEN is based on the algorithm in [9]but with some slight modification. Our experience, in both simulation and real data analysis, suggested that ORIOGEN and the original algorithm in [9]give very similar clustering results. However, the early paper [9]is written in more details and allows our own implementation of the algorithm in Matlab and R. We have carefully examined our coding and found no error.

We believe that the difference between our reported computation time and that in Peddada et al.'s correspondence is mainly due to the efficiency of different computer languages. This is because the bootstrap procedure in Peddada's method requires a very large number of iterations, and JAVA is much more efficient than Matlab or R in looping. A fair comparison of the computational efficiency between two methods should be made on the same platform. Therefore, we implemented both methods in R and reported the computation time in our paper. The computation time of Peddada's method reported on page 3 in [7]is based on our early implementation of the method in Matlab and it used 100,000 bootstrap samples. We also tried ORIOGEN (with 100,000 bootstrap samples), implemented in JAVA by Peddada et al. [8], on our computer to analyze the breast cancer cell line data in [10], and the computation time was 2877 seconds, whereas it only took ORICC (implemented in R) 15.69 seconds. Our computer is a workstation with a 2.30 GHz AMD Athlon(tm) 64 × 2 Dual Core 4400+processor and a 2.00 GB memory. Note that this breast cancer cell line data has been processed and contains just about 1900 genes. Its size is relatively small compared to most current microarray studies. We further applied ORICC and ORIOGEN (with 100,000 bootstrap samples) to a simulated data set that contains 5500 genes. ORIOGEN was implemented using 100,000 bootstraps, with a p-value of 0.0025. The run time for ORICC and ORIOGEN is 74.03 seconds versus 21,104 seconds. We agree with Peddada et al. that the JAVA software ORIOGEN is a very efficient implementation of the algorithm. However, for most microarray studies nowadays involving more than 5,000 genes, without a super powerful computer, it can still take very long for the analyst to obtain the clustering result.

The second comment is that the false positive rates reported in Figure three in [7]were incorrect. We thank Peddada et al. for carefully reading our paper and pointing out this mistake. In our paper, we mistakenly stated the p-value threshold used for Peddada's method. The threshold was 0.5 instead of 0.025. We repeated all simulations in our paper involving Peddada's method using a p-value threshold of 0.025 and the results are presented in Figures ​Figures1,1, ​,2,2, ​,33 and ​and44 at the end of this report. Figures ​Figures1,1, ​,22 and ​and33 are for Simulation 1 in [7]and obtained by imposing the error rate for Peddada's method using threshold 0.025 on Figures ​Figures2,2, ​,33 and ​and44 in [7]. Figure ​Figure44 is for Simulation 2 in [7]and gives Rand's C statistics for the clusters given by Peddada's method using threshold 0.025. As pointed out by Peddada et al. in their correspondence, Peddada's method controls the false positive rate under the nominal level, i.e. the p-value threshold (See Figure ​Figure2).2). However, lower false positive rates are often at the price of increased false negative rates and also higher overall error rates (See Figures ​Figures11 and ​and3).3). Though a threshold of 0.5 seems unreasonable for p-values, it does offer an overall better clustering result than using 0.025 as the threshold. This is further confirmed by Rand's statistics in Figure ​Figure4.4. Except the comparison in false positive rates (Figure ​(Figure2)2) is different from what we stated in our paper, other conclusions in our paper remain unchanged. It is worth noting that Peddada's method can achieve any false positive rate by using the corresponding p-value threshold.

In the situation of clustering microarray data, the clustering result often serves as a hypothesis generating tool, and gives the analysis a more exploratory flavor. Therefore, unlike in the classical hypothesis testing scenario, false negative and false positive are equally important in evaluating the clustering accuracy. Peddada et al. [8,9] treat the clustering problem more in a hypothesis testing way, therefore put more emphasis on controlling the false positive rate. Though our ORICC method is also based on order-restricted inference and has a similar structure to Peddada's method, we view the clustering problem as a model selection problem. By using a consistent model selection criterion, the false positive rate of our ORICC method approaches zero as the number of replicate arrays increases, whereas that of Peddada's method remains around the p-value threshold.

Peddada et al. also had two other comments on our inaccurate description of the optimality of the order-restricted MLE and the test used in Peddada's method. We appreciate their insight and agree that the description should be changed according to their suggestion.

Peddada et al. also reported a new simulation study based on a large number simulated null and non-null gene expressions, in which our one-stage ORICC algorithm was shown to have an inferior performance to the ORIOGEN algorithm. This motivated us to further explore the property of our ORICC algorithms. We repeated the same simulation with a varying number of null genes using ORIOGEN, one-stage ORICC and two-stage ORICC algorithms, and we found that, when null genes consist of the majority of all the genes, two-stage ORICC provides a much more satisfactory performance than one-stage ORICC, and has similar performance to ORIOGEN. See results in Table ​Table2.2. It is worth noting that, in this simulation, the computational advantage of our algorithms remains, especially for the two-stage ORICC. For example, when there are 40000 null genes and 4000 non-null genes, the computational time is 128 minutes, 10.5 minutes and 3.0 minutes for ORIOGEN, one-stage ORICC and two-stage ORICC, respectively, on a workstation with a 2.30 GHz AMD Athlon(tm) 64 × 2 Dual Core 4400+processor and a 2.00 GB memory. Another issue worth mentioning is that ORIOGEN's performance depends on pre-specified p-value cutoff of 0.01 and reclassification p-value of 0.90. On a real data set, finding a proper choice of these cutoffs may be not easy.

To summarize, we thank Dr. Peddada and his colleagues for their insightful discussion, which motivated us to achieve a deeper understanding about the property of ORIOGEN, one-stage ORICC and two-stage ORICC. We hope that our work, their work and discussion, and this report will stimulate further work on clustering microarray data using order-restricted inference.