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?
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=