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.merlevede@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]]

