Dear Wolfgang, dear Simon, Thanks for the explanation, I did not have understand well the purpose of DEXSeq. If you are interested, here are the cases in my dataset where genes have more than 1 exon, but testGeneForDEU did not give a result:
*The simplest case:* res1[1449:1451,] geneID exonID dispersion_CR_est dispersion pvalue EHI_031280:001 EHI_031280 E001 0.25177357 2.517736e-01 0.837492 EHI_031280:002 EHI_031280 E002 0.03580064 3.580064e-02 0.837492 EHI_031280:003 EHI_031280 E003 NA 1.000000e+08 NA padjust EHI_031280:001 0.9746982 EHI_031280:002 0.9746982 EHI_031280:003 NA And the read counts per exon are: exon_counts[1449:1451,] exon_ID FKKM1_1 FPKM1_2 FPKM1_3 FPKM2_1 FPKM2_2 FPKM2_3 1449 EHI_031280:001 9 9 5 33 7 2 1450 EHI_031280:002 155 192 220 453 133 176 1451 EHI_031280:003 0 0 0 0 0 0 One exon in the gene has 0 for all FPKM => 10^8 dispersion and no result for the test. I am not surprised, contrary to *this second case*: res1[1227:1228,] geneID exonID dispersion_CR_est dispersion pvalue padjust EHI_025270:001 EHI_025270 E001 NA 1.00000e+08 NA NA EHI_025270:002 EHI_025270 E002 NA 2.03465e-02 NA NA And the read counts per exon are: exon_counts[1227:1228,] exon_ID FKKM1_1 FPKM1_2 FPKM1_3 FPKM2_1 FPKM2_2 FPKM2_3 1227 EHI_025270:001 0 0 0 0 0 0 1228 EHI_025270:002 11900 3012 12136 10971 7387 5923 Why is there no result for the second exon? Thanks, Jane "Dear Jane, in the case of genes (or better term: loci) with only one exon, DEXSeq's testGeneForDEU does not make sense - it is looking for exons within the locus that behave differently from the average of all the others: there are no others! Hence the parameters are unidentifiable, and the software should not even look at them. We'll robustify the software for this case in future versions (note that the package is still actively being developed). The case of the NA values that are not explained by there only being one exon is more interesting. Can you send us the code to reproduce your problem (this also requires the data, which you can send offline, and in which you can for instance anonymize the gene and sample annotations). Best wishes Wolfgang Hi Jane On 08/08/2011 02:45 PM, Wolfgang Huber wrote: >* in the case of genes (or better term: loci) with only one exon, DEXSeq's *>* testGeneForDEU does not make sense - it is looking for exons within the *>* locus that behave differently from the average of all the others: there *>* are no others! Hence the parameters are unidentifiable, and the software *>* should not even look at them. * Actually, it doesn't. This is why the p value is NA. Also note the column 'testable' in the fData data frame: it is set to FALSE in such cases to indicate that a test was not attempted. You will also find testable=FALSE occasionally in genes with more than one exons, namely whenever all exons or all exons but one have zero counts in all samples. Cheers Simon" 2011/7/29 Jane Merlevede <jane.merlev...@gmail.com> > Hello, > > I am using DEXseq for differential analysis. I have posted some e-mails > about it already on this list, but I have more questions ! > My dataset has 2 conditions with 3 biological replicates per condition: > > pData(ecs) > sizeFactor condition replicate type > Raman_1 0.9638822 nonvir 1 paired-end > HM1_1 1.4000715 vir 1 paired-end > Raman_2 1.0314049 nonvir 2 paired-end > Raman_3 1.0734133 nonvir 3 paired-end > HM1_2 0.7801269 vir 2 paired-end > HM1_3 0.9488409 vir 3 paired-end > > > I found the number of differentially expressed exons. I would like to know > also which are down and which are up regulated. This information is given by > foldChange, which is not provided in DEXSeq (res1). > > res1 <- DEUresultTable(ecs) > head(res1) > geneID exonID dispersion_CR_est dispersion pvalue > EHI_000010:001 EHI_000010 E001 NA 0.02377712 NA > EHI_000130:001 EHI_000130 E001 0.07466945 0.07466945 0.4867354 > EHI_000130:002 EHI_000130 E002 0.12897281 0.12897281 0.9427953 > EHI_000130:003 EHI_000130 E003 0.01960222 0.02031871 0.8432791 > EHI_000130:004 EHI_000130 E004 0.06733720 0.06733720 0.7042977 > EHI_000240:001 EHI_000240 E001 0.04375400 0.04375400 0.7493987 > padjust > EHI_000010:001 NA > EHI_000130:001 0.8639624 > EHI_000130:002 0.9927663 > EHI_000130:003 0.9756197 > EHI_000130:004 0.9503086 > EHI_000240:001 0.9611397 > > So I would like to know, if it could be possible to have baseMean, > baseMeanA, baseMeanB, foldChange, log2FoldCHange, resVarA and resVarB as in > DESeq? > > Since this does not seem to be implemented in DEXSeq, I should start with > the read counts per exon table to get it, but I can't find it neither. I > think it is use in the estimateSizeFactors function, but I don't manage to > access it. It's a pity, I will have to reconstruct the table... > > > Thank you for your help, > Jane Merlevède > [[alternative HTML version deleted]]
_______________________________________________ Bioc-sig-sequencing mailing list Bioc-sig-sequencing@r-project.org https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing