Pairwise comparisons in DEXSeq

[Julien Roux]

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/pipermail/bioconductor/2012-March/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

[areyes]

Hi Julien,

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

Alejandro

[Julien Roux]

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

> 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

[areyes]

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

allExonsPairwise[,-c(1:3)]

[Julien Roux]

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:

> 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 )

[areyes]

Hi!

I meant more like this:

> 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

[Julien Roux]

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?

And I guess your object allExonsPairwise will have only two conditions right?

Correct!

Thanks for your quick answers!
Julien

[areyes]

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