I am analyzing my RNAseq data with two treatments (mock and trt) with three time points each with three replicates using limma. After model fitting etc I calculated qvalues using qvalue package. Then calculated differentially expressed genes for each contrast. They are higher with qvalue cutoff of 0.01 compared p.val cutoff of 0.01 in decideTests (please see the code below).
My major question is how would you incorporate a qval cutoff in your decideTests. I would also like to know why qvalue gives more differentially expressed genes than p.val for the same cutoff since qval should be more stringent.
Appreciate your help.
Sha
#################################################
#My code########
################
#calculating qvals
library(qvalue)
q.val = qvalue(eb.fit$p.val)$qvalues
q.val_contrast1 = q.val [, 'trt_1h - mock_1h']
q.val_contrast2 = q.val [, 'trt_9h - mock_9h']
q.val_contrast3 = q.val [, 'trt_24h - mock_24h']
#Number of DEGs based on qvals
sum(q.val_contrast1 < 0.01)
[1] 6029
sum(q.val_contrast2 < 0.01)
[1] 6008
sum(q.val_contrast3 < 0.01)
[1] 3785
##Differentially expressed genes with decideTests
#####################################
results <- decideTests(eb.fit, p.value = 0.01)
summary (results)
# trt_1h - mock_1h trt_9h - mock_9h trt_24h - mock_24h
# -1 1827 2260 555
# 0 16172 16173 18382
# 1 2932 2498 1994
Something like the code below might work. Assuming that the first coefficient in your
eb.fit
object corresponds to the contrast between treatment and mock at 1 hour:Of course, you can just change
coef
according to the timepoint you're interested in, and you can subsetresults
to remove columns that you don't want to see.