FDR and Pvalue in EdgeR

Question

[saamar.rajput]

I am doing differential expression analysis on my data. I am getting all the results i want except for the FDR value. I want to know both the pvalue for the gene and the adjp for the gene. How to get that ? here is my code and my result

exprs <- as.matrix(read.table(exprsFile, header=TRUE, sep="\t",row.names=1,as.is=TRUE))
libSizes <- as.vector(colSums(exprs))
mobDataGroups <- c("MM", "MM", "MM", "WW", "WW", "WW")
d <- DGEList(counts=exprs,group=factor(mobDataGroups),lib.size=libSizes)
d <- calcNormFactors(d)
d1 <- estimateCommonDisp(d, verbose=T)
d1 <- estimateTagwiseDisp(d1)
fit = glmFit(d1)
result<-glmLRT(fit)

Now the results are

enter code here An object of class "DGELRT"
$coefficients
         (Intercept) y$samples$groupWW
A1BG      -11.676544       -0.24365679
A1BG-AS1  -11.435512        0.29101791
A2M        -6.409277       -0.79089119
A2M-AS1   -14.199847       -0.09916763
A4GALT    -12.529857        1.30770946
14283 more rows ...

$fitted.values
         Donor9.MOCK.2D_S9 Donor10.MOCK.2D_S23 Donor11.MOCK.2D_S37 Donor9.EBOV.2D_S11 Donor10.EBOV.2D_S25 Donor11.EBOV.2D_S39
A1BG            118.953953          102.583303           168.87078          94.028938          106.658704           69.122821
A1BG-AS1        151.409181          130.571986           214.94524         204.387032          231.839861          150.249577
A2M           23086.842583        19909.591121         32774.80868       10561.056642        11979.595154         7763.674009
A2M-AS1           9.428295            8.130756            13.38470           8.600751            9.755986            6.322609
A4GALT           50.606034           43.641543            71.84192         189.107198          214.507672          139.017021
14283 more rows ...

$deviance
     A1BG  A1BG-AS1       A2M   A2M-AS1    A4GALT 
0.7520999 0.6981992 7.4990984 3.6291139 6.1066354 
14283 more elements ...

$method
[1] "oneway"

$unshrunk.coefficients
         (Intercept) y$samples$groupWW
A1BG      -11.677564        -0.2439419
A1BG-AS1  -11.436315         0.2912214
A2M        -6.409282        -0.7908975
A2M-AS1   -14.212585        -0.1006737
A4GALT    -12.532230         1.3094352
14283 more rows ...

$df.residual
[1] 4 4 4 4 4
14283 more elements ...

$design
  (Intercept) y$samples$groupWW
1           1                 0
2           1                 0
3           1                 0
4           1                 1
5           1                 1
6           1                 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$`y$samples$group`
[1] "contr.treatment"


$offset
     [,1]     [,2]    [,3]     [,4]     [,5]     [,6]
x 16.4563 16.30824 16.8067 16.46511 16.59114 16.15739
attr(,"class")
[1] "compressedMatrix"
attr(,"repeat.row")
[1] TRUE
attr(,"repeat.col")
[1] FALSE

$dispersion
[1] 0.09989094 0.07429586 0.06978023 0.42399863 0.16473294
14283 more elements ...

$prior.count
[1] 0.125

$samples
                    group lib.size norm.factors
Donor9.MOCK.2D_S9      MM 14478634    0.9686190
Donor10.MOCK.2D_S23    MM 12627436    0.9577743
Donor11.MOCK.2D_S37    MM 19580908    1.0167715
Donor9.EBOV.2D_S11     WW 14705093    0.9621395
Donor10.EBOV.2D_S25    WW 15048532    1.0664646
Donor11.EBOV.2D_S39    WW 10066760    1.0331801

$prior.df
[1] 10

$AveLogCPM
[1]  2.946438  3.673275 10.224583 -0.359293  3.104285
14283 more elements ...

$table
              logFC    logCPM          LR       PValue
A1BG     -0.3515224  2.946438  0.81165274 0.3676320673
A1BG-AS1  0.4198501  3.673275  1.58309993 0.2083147383
A2M      -1.1410148 10.224583 13.09712089 0.0002957500
A2M-AS1  -0.1430686 -0.359293  0.02828972 0.8664294973
A4GALT    1.8866259  3.104285 13.82666036 0.0002004714
14283 more rows ...

$comparison
[1] "y$samples$groupWW"

$df.test
[1] 1 1 1 1 1
14283 more elements ...

I can see logFC, logCPM, LR and PValue for each gene. I also want to know the adjp for each gene, is that possible?

[Santosh Anand]

you can easily correct Pval according to your favorite method. See ?p.adjust

Answer 1

[Nicolas Rosewick]

extract the p-values from your results and apply FDR correction using p.adjust function :

fdr <- p.adjust(pvalue,"fdr")

Answer 2

[Mark]

Use topTags function to extract the table:

out <- topTags(results, n = "Inf")$table

The default adjustment method is the Benjamini-Hochberg method.