DEXSeq: subsetting ExonCountSet ?
1
0
Entering edit mode
Yvan Wenger ▴ 50
@yvan-wenger-5608
Last seen 6.9 years ago
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
DESeq DEXSeq DESeq DEXSeq • 1.2k views
ADD COMMENT
0
Entering edit mode
@steve-lianoglou-2771
Last seen 22 months ago
United States
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
ADD COMMENT
0
Entering edit mode
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
ADD REPLY

Login before adding your answer.

Traffic: 247 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6