SEQanswers

Go Back   SEQanswers > Sequencing Technologies/Companies > Illumina/Solexa



Similar Threads
Thread Thread Starter Forum Replies Last Post
multiple pairwise comparisons in DESeq and/or edgeR sjm RNA Sequencing 0 03-31-2012 08:12 AM
cuffdiif or cufflinks expr comparisons middlemale Bioinformatics 0 11-05-2010 02:44 AM
R mismatch position in pairwise alignment NicoBxl Bioinformatics 4 10-27-2010 08:36 AM
ChIP-Seq comparisons GERALD Epigenetics 2 07-09-2010 04:36 AM
Performing many pairwise alignments fbarreto Bioinformatics 4 02-27-2010 04:36 AM

Reply
 
Thread Tools
Old 07-12-2012, 07:49 AM   #1
Julien Roux
Member
 
Location: Chicago

Join Date: Dec 2011
Posts: 24
Default Pairwise comparisons in DEXSeq

Dear all,
I am analyzing a RNA-seq experiment including more than 2 conditions. As discussed in this topic, I have created a global exonCountSet, including all conditions and their replicates, to estimate a pooled dispersion from all samples:
https://mailman.stat.ethz.ch/piperma...ch/044094.html

Next, I would like to perform comparisons between specific pairs of conditions included in this exonCountSet object, but I am unsure how to do this... Is this through modified formulas (formula0 and formula1) given to the function testGeneForDEU?

Thanks for your help
Julien
Julien Roux is offline   Reply With Quote
Old 07-13-2012, 08:14 AM   #2
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

Hi Julien,

I think the easiest way is to subset your ecs object and run the pairwise comparisons separately.

Alejandro
areyes is offline   Reply With Quote
Old 07-13-2012, 08:57 AM   #3
Julien Roux
Member
 
Location: Chicago

Join Date: Dec 2011
Posts: 24
Default

Thanks Alejandro,
Yes indeed, that's what I ended up doing:

Code:
> allExonsPairwise <- newExonCountSet(countData = counts(allExons)[, -c(1:3)], design = design(allExons)[-c(1:3),], geneIDs = geneIDs(allExons), exonIDs = exonIDs(allExons))
> allExonsPairwise <- estimateSizeFactors(allExonsPairwise)
> fData(allExonsPairwise)$dispBeforeSharing <- fData(allExons)$dispBeforeSharing
> allExonsPairwise@dispFitCoefs <- allExons@dispFitCoefs
> fData(allExonsPairwise)$dispFitted <- fData(allExons)$dispFitted
> fData(allExonsPairwise)$dispersion <- fData(allExons)$dispersion
> fData(allExonsPairwise)$testable <- fData(allExons)$testable
> allExonsPairwise <- testForDEU( allExonsPairwise, nCores=4 )
Is there a more elegant way to do this?
Julien
Julien Roux is offline   Reply With Quote
Old 07-19-2012, 12:12 AM   #4
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

For the moment I think this us the only option... I will think of a way to be more flexible.

Just one detail, I saw that you are using the dispersion values from the full object in the pairwise comparison... if your objective is to detect pairwise differences, it is better if you recalculate the dispersions just with the two conditions to test. It will increase power because at the moment you are overestimating the variance between the two conditions that you are testing by adding the variance of your other conditions.

Alejandro

ps. you can also subset your dataset just by doing
Code:
allExonsPairwise[,-c(1:3)]
areyes is offline   Reply With Quote
Old 07-19-2012, 07:07 AM   #5
Julien Roux
Member
 
Location: Chicago

Join Date: Dec 2011
Posts: 24
Arrow

Thanks again Alejandro,
If I understood correctly, you suggest to:
- estimate the dispersions on the whole ecs using "estimateDispersions"
- fit the coefficients of the dispersion function from these estimates using "fitDispersionFunction"
- subset the ecs to keep only the conditions you want to compare in a pairwise manner
- reestimate the dispersions on this ecs subset using "estimateDispersions"
- refill fData(ecs)$dispersion with the maximum between the per-gene estimate and the modelled value (second part of "fitDispersionFunction").
Is that correct?

That would translate into this code if I'm correct:
Code:
> allExons <- estimateDispersions( allExons )
> allExons <- fitDispersionFunction( allExons )
> allExonsPairwise <- allExons[, -c(1:3) ]
> allExonsPairwise <- estimateDispersions( allExonsPairwise )
> fData(allExonsPairwise)$dispersion <- pmin(
       pmax( 
          fData(allExonsPairwise)$dispBeforeSharing, 
          fData(allExonsPairwise)$dispFitted,
          na.rm = TRUE ),
          1e8 ) 
> allExonsPairwise <- testForDEU( allExonsPairwise )
Julien Roux is offline   Reply With Quote
Old 07-19-2012, 08:03 AM   #6
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

Hi!

I meant more like this:

Code:
> allExonsPairwise <- allExons[, -c(1:3) ]
> allExonsPairwise <- estimateSizeFactors( allExonsPairWise )
> allExonsPairwise <- estimateDispersions( allExonsPairwise )
> allExonsPairwise <- fitDispersionFunction( allExonsPairwise )
> allExonsPairwise <- testForDEU( allExonsPairwise )
And I guess your object allExonsPairwise will have only two conditions right?

Alejandro
areyes is offline   Reply With Quote
Old 07-19-2012, 09:11 AM   #7
Julien Roux
Member
 
Location: Chicago

Join Date: Dec 2011
Posts: 24
Default

Mmmm, I'm getting confused...
Now I understand you don't agree with Simon and what was discussed here: https://mailman.stat.ethz.ch/piperma...ch/044094.html
(Basically, it is better to estimate the dispersions based on all pooled samples from an experiment, which can include more than the 2 conditions to compare)
Is that correct?

Quote:
Originally Posted by areyes View Post
And I guess your object allExonsPairwise will have only two conditions right?
Correct!

Thanks for your quick answers!
Julien
Julien Roux is offline   Reply With Quote
Old 07-20-2012, 03:29 AM   #8
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

Not necessarily disagree, but I was assuming that you are not expecting the same variability in all your conditions. If this is the case you should maybe estimate them separately, if not, your pooled variance estimation will be more reliable.

I think you could try around and see if they vary much when adding the different conditions, to see how different is the biological variability within each of your conditions.

Alejandro
areyes is offline   Reply With Quote
Reply

Tags
dexseq, exoncountset, pairwise

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:29 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO