Hi, I want to compare DESeq vs DESeq2 and I am getting different number of DEGs which I will expect to be normal. However, when I compare the 149 genes ID that I get with DESeq with the 869 from DESeq2 there are only ~10 genes that are in common which I donât understand (using FDR <0.05 for both). I want to block the Subject effect for which I am including the reduced formula of ~1.
Shouldnât these two methods output similar results? Because at the moment I could interpret my results in different ways. Thanks for your help, Catalina This the DESeq script that I am using: DESeq library(DESeq) co=as.matrix(read.table("2014_04_01_6h_LP.csv",header=T, sep=",", row.names=1)) Subject=c(1,2,3,4,5,1,2,4,5) Treatment=c(rep("co",5),rep("c2",4)) a.con=cbind(Subject,Treatment) cds=newCountDataSet(co,a.con) cds <- estimateSizeFactors( cds) cds <- estimateDispersions(cds,method="pooled-CR", modelFormula=count~Subject+Treatment) #filtering rs = rowSums ( counts ( cds )) theta = 0.2 use = (rs > quantile(rs, probs=theta)) table(use) cdsFilt= cds[ use, ] fit0 <- fitNbinomGLMs (cdsFilt, count~1) fit1 <- fitNbinomGLMs (cdsFilt, count~Treatment) pvals <- nbinomGLMTest (fit1, fit0) padj <- p.adjust( pvals, method="BH" ) padj <- data.frame(padj) row.names(padj)=row.names(cdsFilt) padj_fil <- subset (padj,padj <0.05 ) dim (padj_fil) [1] 149 1 ââââââ library ("DESeq2") countdata=as.matrix(read.table("2014_04_01_6h_LP.csv",header=T, sep=",", row.names=1)) coldata= read.table ("targets.csv", header = T, sep=",",row.names=1) coldata >Subject Treatment >F1 1 co >F2 2 co >F3 3 co >F4 4 co >F5 5 co >H1 1 c2 >H2 2 c2 >H4 4 c2 >H5 5 c2 dds <- DESeqDataSetFromMatrix( countData = countdata, colData = coldata, design = ~ Subject + Treatment) dds dds$Treatment <- relevel (dds$Treatment, "co") dds <- estimateSizeFactors( dds) dds <- estimateDispersions(dds) rs = rowSums ( counts ( dds )) theta = 0.2 use = (rs > quantile(rs, probs=theta)) table(use) ddsFilt= dds[ use, ] dds <- nbinomLRT(ddsFilt, full = design(dds), reduced = ~ 1) resLRT <- results(dds) sum( resLRT$padj < 0.05, na.rm=TRUE ) #[1] 869 [[alternative HTML version deleted]]
______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.