[BioC] DEXSeq: subsetting ExonCountSet ?

[Alejandro Reyes]

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