tweeDEseq glm reduced model specification?

1
@ryan-c-thompson-5618

Last seen 2.1 years ago

Scripps Research, La Jolla, CA

Dear Dr. González,
I am interested in testing tweeDEseq on my own RNA-seq data. I have
noticed that tweeDEseq supports glm fitting for more complex designs
than simple pairwise models. However, it seems that the tweeDEglm
function only supports a null model of "count ~ 1", which makes it
impossible to do do things like test for an interaction term by
comparing "count ~ var1 * var2" against "count ~ var1 + var2". Looking
at the code for tweeDEseq:::anova.glmPT, it looks like the underlying
mechanisms have support fur such comparisons. Are there any future
plans
to support null hypotheses other than "count ~ 1" in the tweeDEglm
function?
Thank you,
-Ryan Thompson
[[alternative HTML version deleted]]

tweeDEseq
tweeDEseq
• 630 views
@gonzalez-ruiz-juan-ramon-6121

Last seen 8.1 years ago

Dear Ryan,
Thank you very much for posting this interesting question.
As you are pointing out, current version of function tweeDEglm only
supports the comparison
count ~ 1
vs
count ~ var1
where var1 encodes for a grouping variable (e.g., cases vs control)
Few changes have to be performed in order to address the comparison
you are proposing
count ~ var1*var2
count ~ var1 + var2
We plan to add this new feature by adding two new arguments in the
function tweeDEglm (modelNull and modelTest).
HOWEVER, you can do such comparison using using glmPT function
(in fact tweeDEglm calls this function gene by gene ans uses lapply to
improve time computation - the idea of encapsulating this repeated
process in a function is also due to the fact that some requirements
have to be full-filled in order to be able to compare both models.
Please, see my comment at the end of this message)
this is an example for analyzing 1 gene (count) for the model you are
interested in:
library(tweeDEseqCountData)
data(pickrell1)
counts <- exprs(pickrell1.eset)
count <- counts[1,]
# simulate two variables for illustrating purposes
n <- ncol(counts)
var1 <- sample(c("A", "B"), n, rep=TRUE)
var2 <- as.factor(sample(c("A", "B"), n, rep=TRUE))
df <- as.factor(data.frame(var1=var1, var2=var2))
mod0 <- glmPT(count ~ var1 + var2, df)
mod1 <- glmPT(count ~ var1 * var2, df)
the p-value corresponding to likelihood ratio test can be obtained by
executing
anova(mod1,mod0)
Therefore, you can use apply or lapply (multicore) for repeating this
process for all genes
The idea of encapsulating this process in tweeDEglm is that by using
'formulas' you can control missing values and guaranty that both
models have the same information (no missing data in counts or var1
var2). If not, the likelihood ratio test is not correct. In addition
we also check whether the models you are comparing are nested. In
other words, you cannot compare
count ~ var1 + var2 +var3
vs
count ~ var1 + var2 + var4
We'll work on that during next days and upload a new version in the
devel 'section' that you will be able to use (if not I can send you
the package including this new feature)
Hope it helps
Best,
Juan
----- Mensaje original -----
De: "Ryan" <rct@thompsonclan.org>
Para: jrgonzalez@creal.cat
CC: "bioconductor" <bioconductor@r-project.org>
Enviados: MiÃ©rcoles, 28 de Agosto 2013 1:20:21
Asunto: tweeDEseq glm reduced model specification?
Dear Dr. GonzÃ¡lez,
I am interested in testing tweeDEseq on my own RNA-seq data. I have
noticed that tweeDEseq supports glm fitting for more complex designs
than simple pairwise models. However, it seems that the tweeDEglm
function only supports a null model of "count ~ 1", which makes it
impossible to do do things like test for an interaction term by
comparing "count ~ var1 * var2" against "count ~ var1 + var2". Looking
at the code for tweeDEseq:::anova.glmPT, it looks like the underlying
mechanisms have support fur such comparisons. Are there any future
plans to support null hypotheses other than "count ~ 1" in the
tweeDEglm function?
Thank you,
-Ryan Thompson
[[alternative HTML version deleted]]

Login before adding your answer.

Traffic: 375 users visited in the last hour