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:
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:
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!!!
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:
Code:
#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), ] )
But, if i change my filtering to:
Code:
resSig <- res[ res$padj < .1, ] up <- subset(resSig, foldChange > 2) down <- subset(resSig, foldChange < 2)
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!!!
Comment