Dear list,
I am not sure the correct way to interpret the nbinomLRT results in
multiple level factors condition.
Here is a toy example. I would to find DE genes after controlling
batch effect in the experiments, in which there are multiple levels.
## make data library(DESeq2)dds <- makeExampleDESeqDataSet(m =
18)colData <- data.frame(row.names = rownames(colData(dds)), sample =
colData(dds)$sample, condition = rep(LETTERS[1:6], each = 3), batch =
factor(rep(c(1, 1, 2), 6)))dds <- DESeqDataSetFromMatrix(counts(dds),
colData = colData, design = ~ batch + condition)dds <-
estimateSizeFactors(dds)dds <- estimateDispersions(dds)## LRTestdds <-
nbinomLRT(dds, reduced = formula(~ batch ))resultsNames(dds)[1]
"Intercept" "batch_2_vs_1" "condition_B_vs_A"
"condition_C_vs_A" "condition_D_vs_A" "condition_E_vs_A"
"condition_F_vs_A"res.1 <- results(dds, name =
"condition_B_vs_A")res.2 <- results(dds, contrast =
list("condition_B_vs_A", "batch_2_vs_1")
Questions:1.) I would like to get the DE genes between B and A, while
controlling for the batch effect. Should I take "res.1" or "res.2", or
both are wrong? and what is the correct way to do it?
2.) Why "res.3 <- results(dds, contrast = c("condition", "A", "b"))"
gave error: "Error in normalizeSingleBracketSubscript(j, x) :
subscript contains invalid names"
3.) In order to get DE genes between conditions not directly listed in
the "resultsNames", is the following codes correct? should
"batch_2_vs_1" be included int the contrast?##for example, DE between
C and F, controlling for batch effect;res.4 <- results(dds, contrast =
list("condition_C_vs_A", "condition_F_vs_A")
Best, Chunxuan
[[alternative HTML version deleted]]

Hi
On 20/08/14 12:43, sh. chunxuan wrote:
> I am not sure the correct way to interpret the nbinomLRT results in
> multiple level factors condition.
Is there a reason why you decided to use 'nbinomLRT' rather than
'nbinomWaldTest', which is suggested in the vignette for standard use
cases?
> Here is a toy example. I would to find DE genes after controlling
> batch effect in the experiments, in which there are multiple levels.
In your case, the LRT allows you to find genes which are affected by
the
experimental condition in _some_ way, i.e., for which you can reject
the
null hypothesis "This gene is expressed at the same level in _all_ six
conditions; i.e., it is not affected by the experimental condition at
all."
However, you want to perform specific contrasts:
> 1.) I would like to get the DE genes between B and A, while
> controlling for the batch effect. Should I take "res.1" or "res.2",
or
> both are wrong? and what is the correct way to do it?
This is much easier with a Wald test. You just run 'nbinomWaldTest'
and
then ask for
results( dds, contrast = c( "condition", "B", "A" ) )
Simon

Hi Chunxuan,
On Aug 20, 2014 12:38 PM, "sh. chunxuan" <hibergo at="" outlook.com="">
wrote:
>
>
>
>
> Dear list,
> I am not sure the correct way to interpret the nbinomLRT results in
multiple level factors condition.
> Here is a toy example. I would to find DE genes after controlling
batch
effect in the experiments, in which there are multiple levels.
> ## make data library(DESeq2)dds <- makeExampleDESeqDataSet(m =
18)colData
<- data.frame(row.names = rownames(colData(dds)), sample =
colData(dds)$sample, condition = rep(LETTERS[1:6], each = 3), batch =
factor(rep(c(1, 1, 2), 6)))dds <- DESeqDataSetFromMatrix(counts(dds),
colData = colData, design = ~ batch + condition)dds <-
estimateSizeFactors(dds)dds <- estimateDispersions(dds)## LRTestdds <-
nbinomLRT(dds, reduced = formula(~ batch ))resultsNames(dds)[1]
"Intercept" "batch_2_vs_1" "condition_B_vs_A"
"condition_C_vs_A"
"condition_D_vs_A" "condition_E_vs_A" "condition_F_vs_A"res.1 <-
results(dds, name = "condition_B_vs_A")res.2 <- results(dds, contrast
=
list("condition_B_vs_A", "batch_2_vs_1")
> Questions:1.) I would like to get the DE genes between B and A,
while
controlling for the batch effect.
You should use the default Wald test for such a comparison, not the
LRT.
See the vignette for an explanation of the difference between the two
tests.
So simply use DESeq() and then use:
results(dds, contrast = c("condition", "B", "A"))
> Should I take "res.1" or "res.2", or both are wrong? and what is the
correct way to do it?
> 2.) Why "res.3 <- results(dds, contrast = c("condition", "A", "b"))"
gave
error: "Error in normalizeSingleBracketSubscript(j, x) : subscript
contains invalid names"
This code gives an error because lowercase "b" is not a level. (It
should
be giving a more understandable error message, not sure why it is not
doing
so here.)
> 3.) In order to get DE genes between conditions not directly listed
in
the "resultsNames", is the following codes correct? should
"batch_2_vs_1"
be included int the contrast?##for example, DE between C and F,
controlling
for batch effect;res.4 <- results(dds, contrast =
list("condition_C_vs_A",
"condition_F_vs_A")
Check the help for ?results and the examples there. You should simply
use:
contrast = c("condition", "C", "F")
Mike
> Best, Chunxuan
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]

Hi Mike, Simon,
Thanks very much for the quick answers.
I moved from DESeq, in which like hood ratio test is used for model
comparison, so the the idea just sticked to my mind.
Later I saw on the other vignette there is a description of the
similar situation specifying "which means that we want to test for the
effect of treatment (the last factor), controlling for the effect of
patients (the first factor)".
Thanks again and best,
Chunxuan
Date: Wed, 20 Aug 2014 14:43:18 +0200
Subject: Re: [BioC] multiple level factors contrast in nbinomLRT,
DESeq2.
From: michaelisaiahlove@gmail.com
To: hibergo at outlook.com
CC: bioconductor at r-project.org
Hi Chunxuan,
On Aug 20, 2014 12:38 PM, "sh. chunxuan" <hibergo at="" outlook.com="">
wrote:
>
>
>
>
> Dear list,
> I am not sure the correct way to interpret the nbinomLRT results in
multiple level factors condition.
> Here is a toy example. I would to find DE genes after controlling
batch effect in the experiments, in which there are multiple levels.
> ## make data library(DESeq2)dds <- makeExampleDESeqDataSet(m =
18)colData <- data.frame(row.names = rownames(colData(dds)), sample =
colData(dds)$sample, condition = rep(LETTERS[1:6], each = 3), batch =
factor(rep(c(1, 1, 2), 6)))dds <- DESeqDataSetFromMatrix(counts(dds),
colData = colData, design = ~ batch + condition)dds <-
estimateSizeFactors(dds)dds <- estimateDispersions(dds)## LRTestdds <-
nbinomLRT(dds, reduced = formula(~ batch ))resultsNames(dds)[1]
"Intercept" "batch_2_vs_1" "condition_B_vs_A"
"condition_C_vs_A" "condition_D_vs_A" "condition_E_vs_A"
"condition_F_vs_A"res.1 <- results(dds, name =
"condition_B_vs_A")res.2 <- results(dds, contrast =
list("condition_B_vs_A", "batch_2_vs_1")
> Questions:1.) I would like to get the DE genes between B and A,
while controlling for the batch effect.
You should use the default Wald test for such a comparison, not the
LRT. See the vignette for an explanation of the difference between the
two tests.
So simply use DESeq() and then use:
results(dds, contrast = c("condition", "B", "A"))
> Should I take "res.1" or "res.2", or both are wrong? and what is the
correct way to do it?
> 2.) Why "res.3 <- results(dds, contrast = c("condition", "A", "b"))"
gave error: "Error in normalizeSingleBracketSubscript(j, x) :
subscript contains invalid names"
This code gives an error because lowercase "b" is not a level. (It
should be giving a more understandable error message, not sure why it
is not doing so here.)
> 3.) In order to get DE genes between conditions not directly listed
in the "resultsNames", is the following codes correct? should
"batch_2_vs_1" be included int the contrast?##for example, DE between
C and F, controlling for batch effect;res.4 <- results(dds, contrast =
list("condition_C_vs_A", "condition_F_vs_A")
Check the help for ?results and the examples there. You should simply
use:
contrast = c("condition", "C", "F")
Mike
> Best, Chunxuan
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]