[BioC] DEXSeq: subsetting ExonCountSet ?

[Steve Lianoglou]

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