Dear Alejandro, Thanks a lot for quick reply! First I appreciate that you pinpoint the typo and will correct it in the analysis. Then I will send you the toBioConductor.RData and the code in a separate e-mail to avoid exceeding the limit of bio-conductor email size. Thanks again, Hui On Jun 26, 2012, at 4:14 AM, Alejandro Reyes wrote: Dear Yao Hui, Thank you for your email and detailed explanations! First, there seems to be a typo in your dispersion formula (which I introduced in the vignette and apparently I forgot to change in the 1.2.0 release). The correct formula would be: formu <- count ~ sample + (type + from) * exon Regarding your question, it sounds like an embarrassing bug... but I could not find in the DEXSeq code a place where the sample names could be switched by mistake. Could you send me your "toBioConductor.RData" object to give it a closer look, as well as the code used to generate the plot?? Best wishes, Alejandro Reyes Dear DEXSeq authors and users, I am using DEXSeq package (1.2.0) to do splicing analysis for human genome. One part of results from the analysis is really confusing me. It seems that for many genes the fitted expression, the fitted splicing and estimated fold changes are inconsistent (actually reversed) with its normalized counts. To clearly explain the problem, I am showing a simple example with only two genes with 30 samples as below. load("toBioConductor.RData") library(DEXSeq) ls() [1] "testgene" ### For this data set, our interest is to investigate the differential usage of exons between two tissue "type", X and N. ### The samples were collected from two batches, So we need to adjust the batch effects in the following model. formu <- count ~ sample + (exon + from)*type testgene <- estimateDispersions(testgene, formula=formu) testgene <- fitDispersionFunction(testgene) f0 <- count~sample + from*exon + type f1 <- count~sample + from*exon + type*I(exon==exonID) testgene <- testForDEU(testgene,formula0=f0, formula1=f1) testgene <- estimatelog2FoldChanges(testgene,fitExpToVar="type" ) res <- DEUresultTable(testgene) n the attached figure, toBioConductor-plotDEXSeq.pdf, for E011, its normalized counts of "N" are clearly smaller than those of "X". However, for both Fitted splicing and Fitted expression, the level of "N" is larger than that of "X". And if we checked the estimated fold change as below, the fold change is consistent with fitted expression and fitted splicing. res["ENSG00000001497:011",] geneID exonID dispersion pvalue padjust ENSG00000001497:011 ENSG00000001497 E011 0.4096852 0.0006160211 0.01334712 meanBase log2fold(X/N) ENSG00000001497:011 5.714405 -1.936449 We also explicitly check the normalized counts for E011 of this gene. As below shows the median of each type. Those of "N" is clearly smaller than "X". dtbl <- data.frame(counts=counts(testgene,normalized=T)["ENSG000000014 97:011",],type=pData(testgene)$type) with(dtbl,tapply(counts,type,median)) X N 8.534784 2.064785 sessionInfo() R version 2.14.0 (2011-10-31) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US LC_NUMERIC=C LC_TIME=en_US [4] LC_COLLATE=en_US LC_MONETARY=en_US LC_MESSAGES=en_US [7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] DEXSeq_1.2.0 Biobase_2.14.0 loaded via a namespace (and not attached): [1] biomaRt_2.10.0 hwriter_1.3 plyr_1.7.1 RCurl_1.91-1 statmod_1.4.14 [6] stringr_0.6 tools_2.14.0 XML_3.9-4 So, have I made any mistakes in the analysis? Many thanks in advance, Hui Hui Yao, Ph.D. Principal Statistical Analyst MD Anderson Cancer Center _______________________________________________ Bioconductor mailing list Bioconductor@r-project.org<mailto:bioconductor@r-project.org> https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor [[alternative HTML version deleted]]