[BioC] How to look at the output of the Poisson data model

[Milena Gongora]

Dear Mark (and edgeR users),

Thanks so much for your feedback. That solved the problem and also was 
good advice regarding the low tag counts.

I am facing another obstacle:
I tried fitting the Poisson model instead of the negative binomial. The 
MAplot and boxplots from this method look good.
Now my problem is, how do I get a table with the resulting 
differentially expressed tags of out of the Poisson model?
topTags() does not work on the object resulting from exactTestNB(), 
because it is not a deDGEList...

My full code is below.
Thanks again!

#-----------------------

d <- new("DGEList", list(data=as.matrix(mydata), group=grouping, 
lib.size=lib_size))
d
head(d$data)
                            A.REP1 A.REP2 B.REP1 B.REP2
100168004_100169791_-       2       4       3       1
100175732_100176961_-       0       3       2       3
100471575_100982094_-       0       0       0       4
101147427_101152326_-       1       0       0       4
101147433_101152326_-       3       1       1       2
101152469_101153215_-       4      21       8      46

qA <- quantileAdjust(d, r.init=1000, n.iter=1)

par(mfrow=c(1,2))
boxplot(as.data.frame(sqrt(d$data)))
boxplot(as.data.frame(sqrt(qA$pseudo)))       

f <- exactTestNB(qA$pseudo, d$group, qA$N * qA$ps$p, r=1000)

#------------------------

Mark Robinson wrote:
> Hi Milena.
>
> Small oversight on my part.  If you use 'as.matrix(mydata)' in place of
> 'mydata', you should have no problems.  This is imposed on the data matrix
> in the next version (which might take a day or so to come online).
>
> A couple of other small points:
>
> 1. You may want to use the new constructor (instead of creating the list
> yourself):
>
> d<-DGEList(data=y,group=grouping,lib.size=lib_size)
>
> 2. In terms of analysis of differential expression, there may not be much
> value in the rows of your table that have VERY low counts (e.g. your
> 100134814_100136947).  In the past, I've noticed some instabilities in the
> NB calculations with very low counts and I've generally filtered out the
> rows with (say) 3 or less total counts.
>
> Cheers,
> Mark
>
>
>
>   
>> Dear edgeR developers and users,
>>
>> Just started using edgeR with next generation sequence (count) data.
>> When calculating alpha, I am getting the following error:
>>
>> alpha <- alpha.approxeb(d)
>> [quantileAdjust] Iteration (dot=1000) 1 :Error in y1/matrix(rep(1,
>> nrow(y1)), ncol = 1) %*% lib.size1 :
>>   non-numeric argument to binary operator
>>
>> I can't work out why I am getting this... any ideas?
>> Thanks!
>>
>> The full code is below:
>> --------------
>> mydata <- read.table("myFile.txt", header=TRUE, row.names="ID")
>> head(mydata)
>>                             A.REP1 A.REP2 B.REP1 B.REP2
>> 100011872_100012727_-       0       0       0       2
>> 100017569_100017878_-       1       0       0       0
>> 100134814_100136947_-       0       0       0       0
>> 100134931_100136947_-       0       2       0       0
>> 100137054_100138100_-       1       0       0       1
>> 100137831_100138100_-       1       0       0       0
>>
>> grouping <- c(1,1,2,2)
>> lib_size <- as.numeric(apply(mydata,2,sum))
>> lib_size
>> [1] 352812 573571 401573 719698
>>
>> d <- new("DGEList", list(data=mydata, group=grouping, lib.size=lib_size))
>> alpha <- alpha.approxeb(d)
>> [quantileAdjust] Iteration (dot=1000) 1 :Error in y1/matrix(rep(1,
>> nrow(y1)), ncol = 1) %*% lib.size1 :
>>   non-numeric argument to binary operator
>> ---------------------
>>
>> --
>> Milena
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor- 
>>     
Milena