DEXSeq log2 fold changes

[yfang01]

I applied the R package DEXSeq to an RNAseq data set and wrote the results to a HTML file. Comparing the expression plot and log2 fold change in the corresponding result table, the log2foldchange values on the table seem conservative relative to those shown on the plot. Could anyone help me to understand this discrepancy?

For example, in an analysis of the pasilla sample data set, exon 10 in FBgn0010909 looks to be ~4 fold lower in the treated data:

http://www.cgr.liv.ac.uk:8088/illum/ID0095DEU_2f847c24e9311da6/model1/files/FBgn0010909expression.html

but in the results table it has a log2foldchange value of 0.64:

http://www.cgr.liv.ac.uk:8088/illum/ID0095DEU_2f847c24e9311da6/model1/files/FBgn0010909results.html

The following commands were put in to do the analysis:

library("DEXSeq")
library("pasilla")
data("pasillaExons", package="pasilla")
pasillaExons <- estimateSizeFactors(pasillaExons)
pasillaExons <- estimateDispersions(pasillaExons)
pasillaExons <- fitDispersionFunction(pasillaExons)
pasillaExons <- testForDEU( pasillaExons,nCore=8)
pasillaExons <- estimatelog2FoldChanges( pasillaExons )

ecs1 <- pasillaExons
#
formuladispersion <- count ~ sample + ( condition + type ) * exon
pasillaExons <- estimateDispersions( pasillaExons, formula = formuladispersion )
pasillaExons <- fitDispersionFunction(pasillaExons)
formula0 <- count ~ sample + type * exon + condition
formula1 <- count ~ sample + type * exon + condition * I(exon == exonID)
pasillaExons <- testForDEU( pasillaExons, formula0=formula0, formula1=formula1,nCore=8)
pasillaExons <- estimatelog2FoldChanges( pasillaExons )
ecs2 <- pasillaExons

Then the results were written by:
DEXSeqHTML( ecs1, FDR=0.1, color=c("#FF000080", "#0000FF80"),path="/agf/illum/ID0095DEU_2f847c24e9311da6/model1/",file="pasilla_DEU_res.html")

#
DEXSeqHTML( ecs2, FDR=0.1, color=c("#FF000080", "#0000FF80"),path="/agf/illum/ID0095DEU_2f847c24e9311da6/model2",file="pasilla_DEU_res.html")

The results are available from:

http://www.cgr.liv.ac.uk:8088/illum/ID0095DEU_2f847c24e9311da6/

Thanks in advance if anyone can give a reasonable explanation.

[areyes]

Dear yfang01,

Thanks for your question and interest in DEXSeq. The reason is that the fold change that DEXSeq gives is not calculated directly from the counts, what DEXSeq does is:

1. The interaction coefficient is taken from the model: count ~ condition * exon. (this is in order to distinguish changes in gene expression levels from changes in the usage of individual exon, imagine a situation in which a whole gene is differentiallty expressed, your fold changes would mainly reflect changes in expression, not in exon usage):

coefficients <- DEXSeq:::fitAndArrangeCoefs(pasillaExons, "FBgn0010909", frm=count~condition*exon )
effects <- DEXSeq:::getEffectsForPlotting(coefficients, "condition", averageOutExpression=TRUE)

2. The last variable "effects" will contain the natural logarithm of the "exon usage" effects, which should approximate the real counts. We next apply the variance stabilising transformation to this coefficients (described in the DEXSeq paper) and then transform them to the log2 scale.

effects <- log2( DEXSeq:::vst( exp( t( effects ) ), pasillaExons ) )
foldChange <- effects[,"untreated"] - effects[,"treated"]
> foldChange
E001 E002 E003 E004 E005 E006
-0.002012139 -0.077514494 -0.085524653 -0.086438377 -0.231731049 -0.168603093
E007 E008 E009 E010 E011 E012
-0.043928621 0.014562948 0.036142063 0.641266891 0.017428325 0.032003742
E013 E014 E015 E016 E017 E018
-0.191428463 0.003098754 0.036862405 0.021424340 -0.062929341 -0.048894328
E019 E020 E021 E022 E023
-0.041459857 -0.088129731 -0.012760163 -0.093161888 -0.135868986


For check the advantages of this fold change have a look at the DESeq vignette!

Alejandro