False discovery rate from permutation testing?

Question

I have conducted a search for genetic interactions using a simple dosage model:

Y ~ A + B + AB

where Y is the phenotype, in this case, gene expression values and A and B are vectors of genotype information for ~500 samples. I wish to determine a signficance threshold using permutation testing in order to correct for multiple testing.

To date, I have recalculated the p-values for the interaction term (AB) for 100 permutations (I permuted the phenotype values) and am unsure how to proceed in order to derive a false discovery rate (FDR).

Any suggestions?

Thanks, D.

Answer 1

What I know that people do (and I do it myself in some sort of way in GWAS studies) is that you combine all your permuted p values into a null distribution and then just see how many p values are above a threshold in your real experiment and your permuted null.

So your FDR would then be something like: number of null p values < x / number of real p values < x

Answer 2

1) I agree with suncoolsu that FDR is not estimated (to the best of my knowledge, which is slim), but controlled for.

2) Once you have p values you can use something like - p.adjust(my_p_values, method = "BH"). And you will get the adjusted p values using the Benjamini Hochberg procedure. Any p value that is bellow (let's say) 0.05, can be rejected for a Q=0.05.

Answer 3

The adjusted analogue to the p-values you are probably looking for is the q-value, which is described by Storey as:

(...) [giving] the scientist a hypothesis testing error measure for each observed statistic with respect to pFDR.

The p-value accomplishes the same goal with respect to the type I error, and the adjusted p-value with respect to FWER.

So, q-value is to FDR as adjusted p-value is to FWER.

In R, there is a qvalue package which can produce these estimates given p-values.

Answer 4

An old question but anyway: I've just been looking into this since it seems to be one of these many things that people do in genomics without explicitly saying what they are doing or why it's justified.

The idea is laid out pretty clearly in this paper: http://bioinformatics.oxfordjournals.org/content/21/23/4280.full where they also point out that you really need an estimate of $\pi_0$, the proportion of negatives. Estimating $\pi_0$ is often done using the Storey method in the qvalue R package, but if you're doing permutations that implies you don't think your p-values are well-calibrated (Storey assumes calibrated p-values). Then the only really safe procedure is to assume $\pi_0=1$ (what Sander suggests implicitly above), which is at least conservative.

As an aside, there is no valid permutation strategy for testing interaction terms.