R: How to change the default FDR value imbedded in the nbinomTest() function?

688 Views Asked by At

I have a table, here's the start:

WE1_counts  WE2_counts  M1_counts   M2_counts   M3_counts
YAL008W 465 291 911 926 946
YBR255W 1040    1357    1428    1304    1112
YGR131W 95  170 230 113 138
YNL003ßC    1800    3107    3979    3012    2899
YBR135W 1094    2143    936 860 561

The goal is to find differentially expressed genes between conditions "W" and "M". Here's my code to do so:

source("http://bioconductor.org/biocLite.R")
biocLite("DESeq")

library("locfit")
library("lattice")
library("DESeq") 

condition= c("W", "W", "M", "M", "M")
#data = table given above
cds1 = newCountDataSet(data, condition)

#Normalization
cds1=estimateSizeFactors(cds1)
sizeFactors(cds1)
head( counts( cds1, normalized=TRUE ) )

#Estimate dispersion
cds1 = estimateDispersions( cds1 )
dev.new()
plotDispEsts( cds1 )
head( fData(cds1) )

# Plotting mean vs log2FoldChange
res = nbinomTest( cds1, "W", "M")
dev.new()
plotMA(res)

My question is this: the nbinomTest() function gives me back differentially expressed genes based on a 10% FDR. Is there way to change this number? Can I check for differentially expressed genes at, for example, 5% FDR?

1

There are 1 best solutions below

0
On

nbinomTest merely adjusts for multiple comparisons, you then apply the cutoff yourself.

counts(cds1)[which(res$padj < my.fdr),]

If on the other hand, you wan't to control the cutoff (specifically, the α \alpha) of the FDR procedure itself, then you might need this (have'nt used it tho)

PS: BTW, while coding in R, please use consider using <- instead of =