Significant up-regulated and downregulated genes with DESeq

[alisrpp]

Hi,

I'm trying to identify the significant up-regulated and down-regulated genes in 6 different moments of the life cycle of "my" animal.

I have 11 libraries (2 replicates/stage, except for one stage that i couldn't sequence a replicate). What i'm doing is to compare one stage with the stage immediately after using DESeq.

This is my script:

#Load data
Cele_SPvsEB = read.csv (file.choose(), header=TRUE, row.names=1)
CeleDesign <- data.frame(
  row.names = colnames(Cele_SPvsEB),
  condition = c("SP", "SP", "EB", "EB"))

#to load conditions w/out pData file: conditions
conds = factor(c("SP", "SP", "EB", "EB"))

#Make cds file
cds <-newCountDataSet(Cele_SPvsEB, conds )

#Est size factor = normalize for library size
cds<-estimateSizeFactors(cds) 

#estimate dispersion from the expected
cds <- estimateDispersions(cds, method=c("pooled"), fitType=c("local"))
str(fitInfo(cds))

#plot dispersion
plotDispEsts <- function( cds )
{
   plot(
      rowMeans( counts( cds, normalized=TRUE ) ),
      fitInfo(cds)$perGeneDispEsts,
      pch = '.', log="xy", ylab="dispersion", xlab="mean of normalized counts" )
   xg <- 10^seq( -.5, 5, length.out=300 )
   lines( xg, fitInfo(cds)$dispFun( xg ), col="red" )
}
plotDispEsts( cds )

#call differential expression
res <- nbinomTest( cds, "SP", "EB" )

#there is no expression for a gene accross all samples... need to remove those genes
whichIsNA<-which(is.na(res[,8]))
res<-res[-whichIsNA,]

#plot volcano
plotDE <- function( res )
   plot(
      res$baseMean,
      res$log2FoldChange,
      log="x", pch=20, cex=.3,
      col = ifelse( res$pval < .05, "red", "black" ) )

plotDE( res )

#filter for upregulated and downregulated genes

resSig <- res[ res$pval < .05, ]
up <- subset(resSig, foldChange > 2)
down <- subset(resSig, foldChange < 2)

#list most sig genes
head( resSig[ order(res$pval), ] )

By using this script i found around 600 genes in both up-regulated and down-regulated categories for this example (i.e. SP vs EB).

But, if i change my filtering to:

resSig <- res[ res$padj < .1, ]
up <- subset(resSig, foldChange > 2)
down <- subset(resSig, foldChange < 2)

Then i have 0 genes up or down-regulated.

If i just keep going only filtering by pval < 0.05, is that enough to consider my genes significantly differentially expressed?

Do i have any error in my script?

Any suggestions?

Thanks!!!

[dpryan]

If i just keep going only filtering by pval < 0.05, is that enough to consider my genes significantly differentially expressed?

No. No one will (or at least should) believe those are differentially expressed without additional evidence. Also, keep in mind also that that's probably log2(fold-change), so your up/down regulated may not mean what you think (I haven't used DESeq1 in a while, but I remember this being the case). So, things to think about trying are as follows:

1. Try DESeq2, which has a number of very nice improvements.
2. Use the genefilter package to help properly perform independent filtering. This can drastically up your power.

These two options aren't mutually exclusive, of course.

[sindrle]

Can you calculate fold change from log2FC and the filter out above say 1.5?

You set it at 2, why not lower? cant you set it at 1.5?

[dpryan]

Can you calculate fold change from log2FC and the filter out above say 1.5?

Yes, that works.

You set it at 2, why not lower? cant you set it at 1.5?

The thresholds are all semi-arbitrary, you could set them to whatever makes sense for your experiment.

[sindrle]

How do I know "it makes sense". Say you have a sign. diff. exp. gene at about 30 FPKM and a FC of 1.5. Maybe not much? But its biological function makes a lot of sense in the context.
What then do you do?

[dpryan]

There's no rule I can give you for determining "what makes sense", it'll depend entirely on the system you're using, your end goals, the resources (financial and otherwise) at hand, and how risk averse you are. At some point you just have to become comfortable with contextualizing results in terms of biological implications (I only filter by adjusted p-value for this reason). Back in the early microarray days, using fold-change thresholds was pretty common, seemingly due to the methods not being that great. These days, I haven't seen as many people doing that.

[sindrle]

Yeah, that was what i thought!

Question about "filtering", does it make sense only to keep the intersect of Cuffdiff, edger and deseq2? Or how to know which results to use?

[dpryan]

Presumably you're going to do some validations via qPCR or something. So just include a few unique hits from each in an initial screen to get an idea which is better fitting your experiment. I would expect either edgeR or DESeq2 to produce the best results, but it's easy enough to check.

[sindrle]

Yes, I have calculated the result of one gene so far with qPCR. And edgeR is in the lead!

Make that 3/3 so far.