Validity of replacing q-value with p-values from multiple testing of different tools

In this Science paper, the authors are partly interested in trying to identify genes that cycle (that oscillate through the day in terms with the biological clock). Here's an extract from the methods explaining their method:

Time Series Analysis for Circadian Cycling:

RNA cycling was assessed by three programs, COSOPT (35), JTK cycle (75) and ARSER (76) with expression level cutoffs, > 0.05 RPKM for intron, > 0.5 RPKM for exon and > 0.0625 RPKM for antisense. For COSOPT and JTK cycle analyses, data was detrended by linear regression. A cycling gene was considered if two out of three programs detected cycling with threshold of p<0.05. The period and phase from ARSER were used for further analysis. ForChIP-seq peak analysis, two cycles were concatenated and the cycling was analyzed with ARSER (p<0.05).

My question is, how much is it statistically acceptable/valid to do multiple testing (over 20000 tests) using different softwares and then choosing those that are significant at p-value (not FDR) < 0.05 from all the softwares?

Hard to say, since different methods can give results that are more or less correlated.

Assuming that there are no cycling genes (H0), and that the predictions and genes are independent (both assumptions being probably not true), we can calculate the expected number of false positives. For a given gene, the probability of a FP in a test is 0.05. The probability that at least two out of three give a hit is then $0.05^3 + 3\cdot0.05^2\cdot0.95 = 0.00725$. In a set of 20,000 genes I would thus expect 145 false positives at the given $p$-value thresholds.

This is much less then what the authors detected (if I'm not mistaken). However, were the methods perfectly correlated, the expected number of false positives would be around a thousand. Are you sure there were only 20,000 tests, though? Without going into the details of the paper, they look at both, different exons and different introns, and 20,000 is just the rough number of the genes, not accounting for exons and introns.

• January, thanks for your answer. I chose this 20K genes arbitrarily (as an example of the amount of testing that happens). I now see that this method may not be as bad as I imagined, but could you explain a bit how you get this formula (that gives 0.00725)?
– Arun
Jul 4, 2013 at 14:46
• Binomial distribution with $n=3, k=2$ and $p= 0.05$ :-) $0.05^3$ - probability that all tests are significant; $0.05 \cdot 0.05 \cdot 0.95$ -- probability that first two are significant etc. Jul 4, 2013 at 16:38
• Thank you! To summarise, the probability that a FP is observed at 5% significance with the intersection of at least 2 out of 3 programs is 0.00725 (assuming independence). This allows me to calculate it for the actual data. I was wondering how one could ascertain the validity of the test. Your answer really helped that point of mine I was interested in. After looking at the answer, it looks so obvious to me! Just don't seem to come up with it myself! :)
– Arun
Jul 4, 2013 at 19:07

The conservative answer is clearly no- you should not.

To take this to an extreme- assuming you have infinitely many different software, each provides type I error control over the 20,000 genes, but any two software differ slightly. You will then always reject your global null hypothesis even if it is true.

Multiple software is indeed a multiple testing problem. The good news is that as @January pointed out- since you would expect strong correlations over different software, the type $I$ error rate inflation is slow, with the growing number of software. The power differences between the software might however be large, so tying out a small number of different software might give you a large power gain at the cost of a small type $I$ error inflation so that it could actually be justified.

• The problem is that the reason they go for different softwares is that there's normally a huge non-overlap between these softwares (meaning the correlation will be typically low). This is also my experience with these softwares on my data. But I choose one and fetch results conservatively. I am beginning to get the feeling that, in those cases, from what you and January point out, it'll most definitely fail...?
– Arun
Jul 5, 2013 at 9:11
• Assuming all software provide type $I$ error control over the 20,000 genes, but in very different manners, so that their global p-value are uncorrelated (under the global null), you will have fast type $I$ error inflation. I believe however, this is implausible as each field has a small number of favorite error control methods. You could use simulated data to quantify the agreement over different software under the null. Jul 5, 2013 at 9:15