Entering edit mode
7.1 years ago
saamar.rajput
▴
80
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?
you can easily correct Pval according to your favorite method. See
?p.adjust