Fdr Q Values Calculations Using Q-Value Package

Question

[viavialin]

Hi Everyone,

I'm trying to determine the q values of my dataset using p-values generated by Fisher's tests from PLINK (allelic and genotypic). The p values are associated with SNPs relevant against clinical outcome. I am trying to determine the most significant SNPs from my analysis out of 70 SNPs.

After comparing both 'bootstrap' and 'smoother' qvalues output. I realize that some of the qvalues are < original p values after a certain threshold which is not the value reporting "qobj$pi0".

This output is definitely odd since I should expect qvalues to be > p values, especially when the pvalue is 1 and the qvalue (bootstrap) reported is 0.17286 (qobj$pi0 value)

If anyone can give me some clues, that would be great! I would be happy to paste my p-values if necessary.

Thanks in advance.

POST EDITED to include p-values (not sure how to attach a file here)

0.006211 0.231 0.1929 0.2727 0.0128 0.2003 0.02752 0.006144 1 0.6376 1 0.003863 0.002481 0.05531 0.02465 0.0141 0.0721 0.01738 0.0713 0.01848 0.02205 0.1721 0.03096 0.02307 0.06568 0.004345 0.0006609 0.00259 0.005319 0.0004572 0.0143 0.001896 0.02808 0.01518 0.1623 0.184 0.2339 0.2514 0.07164 0.04881 0.09818 0.1057 0.08894 0.07172 0.09818 0.1057 0.09542 0.1036 0.2362 0.2519 0.2516 0.2512 0.2186 0.1833 0.2036 0.1829 0.07978 0.04875 0.09542 0.1036 0.09361 0.1039 0.3349 0.3466 0.093 0.1488 0.042 0.02281 0.1044 0.07099 0.8093 0.8393 0.4839 0.441 0.8658 1 0.2296 0.2512 0.1019 0.1296 0.1154 0.1023 0.3635 0.1534 0.2095 0.2528 0.006035 0.001488 0.01061

[David Quigley]

Post the P values. Note that qvalue does not work as expected unless your P values have the correct distribution: uniform, with a possible (hopefully present) over-representation close to zero.

Answer 1

[Charles Warden]

1) I think you are using the wrong result: qobj$qvalues is the array of q-values. pi0 is a parameter used to estimate the q-values (explained in more detail in the q-value paper: http://www.pnas.org/content/100/16/9440.short)

2) For SNP studies, I think a Bonferroni correction might be more common than an FDR correction (at least that is what is used by the GWAS catalog). This is also obviously much easier to calculate, but should be more stringent.

[viavialin]

Thanks for your insights. From what I understood in that paper, pi0 is an estimate of proportion of insignificant tests in your dataset, and that the q value is a version of FDR corrected p-values. Perhaps I am understanding this inaccurately.

I agree that Bonferroni correction is the more commonly used tool but since I only have 70 SNPs to query, I am already getting less than a handful of positive SNPs.

Could it be that I am using qvalue incorrectly? The command I am using in the qvalue package, qobj <-qvalue(p,pi0.meth="bootstrap")

On the other hand, I tried p.adjust option in R with p.adjust(p,method="fdr") and obviously it uses the Benjamini and Hochberg to calculate the q values. The corrected q values appear to be fine. Do you suggest that I stick to the B&Y method instead of the Storey method?

Thanks for your patience with my questions, I am fairly new to statistics and genetic studies.

[Charles Warden]

I think you are using q-value correctly, you just need to recognize that the q-values are the end result (not the pi0 value). I also don't think you need to specify any special parameters.

For example, see the demo code on the qvalue website:

http://www.bioconductor.org/packages/release/bioc/vignettes/qvalue/inst/doc/qvalue.R

To make my life easier, I typically do use p.adjust. However, the Storey q-value is often a little lower than the B-H q-value. So, sometimes I do find reporting q-values to be easier.

Perhaps more importantly, how is it that you only have 70 SNPs? All the SNP arrays that I know of have way more (and I would expect NGS data to also typically have way more). By "way more", I mean somewhere roughly around 10^6 probes (which is why you would often use the Bonferroni correction).