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
<catagui@gmail.com>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]]
>
>
> ______________________________________________
> 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.
>
>
[[alternative HTML version deleted]]