You're on the wrong list. This is more appropriate on the bioconductor mailing list.
On Mon, May 5, 2014 at 9:42 AM, Catalina Aguilar Hurtado <[email protected]>wrote: > 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]] > > > ______________________________________________ > [email protected] 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. > > [[alternative HTML version deleted]]
______________________________________________ [email protected] 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.
