FDR and Pvalue in EdgeR
2
0
Entering edit mode
7.1 years ago

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?

RNA-Seq rna-seq R • 10k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
7.1 years ago

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

fdr <- p.adjust(pvalue,"fdr")
ADD COMMENT
2
Entering edit mode
6.9 years ago
Mark ★ 1.6k

Use topTags function to extract the table:

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

The default adjustment method is the Benjamini-Hochberg method.

ADD COMMENT

Login before adding your answer.

Traffic: 1541 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