How to correctly calculate Storey-Tibshriani FDR Q Values using R Package "QValue"?
0
0
Entering edit mode
7.9 years ago
Tom ▴ 50

Hi all. I am trying to replicate the results of another paper.

The aim: I have a P value. Then I randomised the data set 1,000 times, and now I have another set of randomised P Values. As an example, let's say I have the p value below, and I've randomised the data set 10 times and I have 10 randomised p values:

Actual P Value = 0.03
Random P Values:

0.06
0.21
0.003
0.21
0.017
0.21
0.48
0.21
0.48
0.21

The exact technique that I want to do is: "Calculate the false discovery rate (FDR) (Q), defined as the number of expected false positives over the number of significant results (Storey and Tibshirani, 2003)".

I wrote this code using python, rpy2 and R:

import sys
import rpy2
import rpy2.robjects as ro
from rpy2.robjects import r
from rpy2.robjects.packages import importr

devtools = importr("devtools")
qvalue = importr("qvalue")

list1 = [float(line.strip()) for line in open(sys.argv[1])] #This is a list of 1,000 floats (randomised P Values)
pvals = ro.FloatVector(list1)
rcode = 'qobj <-qvalue(p=%s)' %(pvals.r_repr())
res = ro.r(rcode)
r_output1 = 'qobj$pvalue'
r_output2 = 'qobj$qvalue'
r_pvalue = ro.r(r_output1)
r_qvalue = ro.r(r_output2)
ListOfValues = zip(r_pvalue,r_qvalue)
for i in ListOfValues:
            print str(i[0]) + "\t" + str(i[1])

The output looks like this (a snapshot of my output for 1,000 p values):

P value Q value
0.21    0.032
0.48    0.04
0.21    0.032
0.81    0.068
0.017   0.017
0.48    0.047
0.81    0.06
0.017   0.017
0.211   0.03
0.81    0.068
0.069   0.023
0.069   0.023
0.069   0.023
0.48    0.047

I do not want an individual Q value for every randomised P value. I have an original P value (in this case 0.03), and I want one Q value to be associated with this P value; based on the set of randomised P Values. Could someone tell me where I've gone wrong in this code; or how I extract this information?

Thanks.

r fdr python q value p value • 4.1k views
ADD COMMENT
0
Entering edit mode

Hi, I am not an expert, but I have the feeling that you mix fdr and bootstrap. May be those articles could help you with FDR. http://www.nature.com/nmeth/journal/v11/n4/full/nmeth.2900.html http://www.nature.com/nbt/journal/v27/n12/full/nbt1209-1135.html Best.

ADD REPLY

Login before adding your answer.

Traffic: 1800 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6