Hello everybody,
I have 11 libraries representing 6 different conditions (some are in
triplicates). Is there a simple way to test for differential
expression using DEXSeq between two libraries I?actually choose?
For example, in the DESeq pacakge, we use the command res =
nbinomTest( cds, "untreated", "treated" ) which clearly specifies the
conditions between the two samples to test, but I do not see any
possibility to input the conditions I want in the DEXSeq function
testForDEU(). I would like to be able to use dispersions calculated
with all my conditions even if I?test only condition pairs (like in
DESeq basically).
Many thanks for your help,
Yvan
Hi Yvan,
On Tue, May 7, 2013 at 1:26 AM, Yvan Wenger <yvan.wenger at="" unige.ch="">
wrote:
> Hello everybody,
>
> I have 11 libraries representing 6 different conditions (some are in
> triplicates). Is there a simple way to test for differential
> expression using DEXSeq between two libraries I actually choose?
>
> For example, in the DESeq pacakge, we use the command res =
> nbinomTest( cds, "untreated", "treated" ) which clearly specifies
the
> conditions between the two samples to test, but I do not see any
> possibility to input the conditions I want in the DEXSeq function
> testForDEU(). I would like to be able to use dispersions calculated
> with all my conditions even if I test only condition pairs (like in
> DESeq basically).
This is relatively simple to do. You first have to subset out the
samples that belong to your condition of interest then run the
`testForDEU`. You also have to be careful to remove the levels of your
"condition" factor that have been dropped due to your subsetting your
ExonCountSet object.
Let's use the `pasillaExons` ExonCountSet from the pasilla package to
see how to test only the samples of `type == "single-read"`
R> library(pasilla)
R> data("pasillaExons")
R> design(pasillaExons)
condition type
treated1fb | treated single-read
treated2fb | treated paired-end
treated3fb | treated paired-end
untreated1fb| untreated single-read
untreated2fb| untreated single-read
untreated3fb| untreated paired-end
untreated4fb| untreated paired-end
## Now take only samples of `type == "single-read"`
R> some <- pasillaExons[,design(pasillaExons)$type == "single-read"]
R> design(some)
condition type
treated1fb | treated single-read
untreated1fb| untreated single-read
untreated2fb| untreated single-read
## The problem is that the `paired-end` read level is still
## in the `type` column of the design:
R> design(some)$type
[1] single-read single-read single-read
Levels: paired-end single-read
## This will trip you up when you run the `testForDEU` as the
## model matrix will have 0-columns that are generated from
## the levels that have been removed from the design. This
## is easy to fix:
R> design(some) <- droplevels(design(some))
## Now let it rip
R> result <- testForDEU(some, ...)
Alejandro: perhaps adding a `"[","ExonCountSet"` method to do this
automatically would be useful, eg. in the R/class_and_accessors.R file
you could have something like:
setMethod("[", "ApaCountSet", function(x, i, j, ..., drop = FALSE) {
x <- callNextMethod()
design(x) <- droplevels(design(x))
x
})
HTH,
-steve
--
Steve Lianoglou
Computational Biologist
Department of Bioinformatics and Computational Biology
Genentech
Dear Steve,
Thanks for the suggestion!
I have included this function into DEXSeq.
Alejandro
> Hi Yvan,
>
> On Tue, May 7, 2013 at 1:26 AM, Yvan Wenger <yvan.wenger at="" unige.ch=""> wrote:
>> Hello everybody,
>>
>> I have 11 libraries representing 6 different conditions (some are
in
>> triplicates). Is there a simple way to test for differential
>> expression using DEXSeq between two libraries I actually choose?
>>
>> For example, in the DESeq pacakge, we use the command res =
>> nbinomTest( cds, "untreated", "treated" ) which clearly specifies
the
>> conditions between the two samples to test, but I do not see any
>> possibility to input the conditions I want in the DEXSeq function
>> testForDEU(). I would like to be able to use dispersions calculated
>> with all my conditions even if I test only condition pairs (like in
>> DESeq basically).
> This is relatively simple to do. You first have to subset out the
> samples that belong to your condition of interest then run the
> `testForDEU`. You also have to be careful to remove the levels of
your
> "condition" factor that have been dropped due to your subsetting
your
> ExonCountSet object.
>
> Let's use the `pasillaExons` ExonCountSet from the pasilla package
to
> see how to test only the samples of `type == "single-read"`
>
> R> library(pasilla)
> R> data("pasillaExons")
> R> design(pasillaExons)
> condition type
> treated1fb | treated single-read
> treated2fb | treated paired-end
> treated3fb | treated paired-end
> untreated1fb| untreated single-read
> untreated2fb| untreated single-read
> untreated3fb| untreated paired-end
> untreated4fb| untreated paired-end
>
> ## Now take only samples of `type == "single-read"`
> R> some <- pasillaExons[,design(pasillaExons)$type == "single-read"]
> R> design(some)
> condition type
> treated1fb | treated single-read
> untreated1fb| untreated single-read
> untreated2fb| untreated single-read
>
> ## The problem is that the `paired-end` read level is still
> ## in the `type` column of the design:
> R> design(some)$type
> [1] single-read single-read single-read
> Levels: paired-end single-read
>
> ## This will trip you up when you run the `testForDEU` as the
> ## model matrix will have 0-columns that are generated from
> ## the levels that have been removed from the design. This
> ## is easy to fix:
> R> design(some) <- droplevels(design(some))
>
> ## Now let it rip
> R> result <- testForDEU(some, ...)
>
> Alejandro: perhaps adding a `"[","ExonCountSet"` method to do this
> automatically would be useful, eg. in the R/class_and_accessors.R
file
> you could have something like:
>
> setMethod("[", "ApaCountSet", function(x, i, j, ..., drop = FALSE) {
> x <- callNextMethod()
> design(x) <- droplevels(design(x))
> x
> })
>
> HTH,
> -steve
>
> --
> Steve Lianoglou
> Computational Biologist
> Department of Bioinformatics and Computational Biology
> Genentech
