[BioC] DEXSeq on two-exon genes: how to specify a formula without redundant terms

Narayanan, Manikandan (NIH/NIAID) [E] manikandan.narayanan at nih.gov
Fri May 24 17:08:35 CEST 2013


Great, thank you Alejandro - I will check out the code and let you know how it goes. 

In the mean time, I also found a small typo in the DEXSeqHTML function. The line 
  results[which(is.na(results$pvalue)), ][, c("pvalue", "padjust")] = 1
breaks with an error if there are no NA pvalues, and can be patched by changing to:
    results[which(is.na(results$pvalue)),c("pvalue", "padjust")] = 1 

While you are at it, it would also be useful if extraCols is indexed by geneIDs as rownames - that way we can send in extra info for as many 
genes in any order and it would be properly added to genetable for display in the DEXSeqHTML report. In detail, if you could change:
  genetable <- cbind(extraCols, genetable)
to 
  genetable <- cbind(extraCols[match(genetable$geneID, rownames(extraCols)),], genetable]) 
or something similar and change the documentation appropriatley, that would be quite helpful to users too!


Thanks, 
Mani



________________________________________
From: Alejandro Reyes [alejandro.reyes at embl.de]
Sent: Friday, May 24, 2013 8:19 AM
To: Narayanan, Manikandan (NIH/NIAID) [E]
Cc: bioconductor at r-project.org
Subject: Re: [BioC] DEXSeq on two-exon genes: how to specify a formula without redundant terms

Dear Manikandan Narayanan,

Sorry for the delayed reply!

Thanks for your e-mail! It has led us to spot two minor bugs and
therefore modified the DEXSeq code in two aspects:

Firstly, redundant terms were removed in DEXSeq in a very naive way:  a
lm.fit model was first fitted and the columns from the model matrix with
NA coefficients were removed from the model matrix, now a proper
function has been written for this purpose (with basically the code to
what you also proposed). Secondly, our "if" statement that you mentioned
was wrongly located within the function, which was causing your error
message. It has been relocated and now it seems to work for its purpose
( I tested with your model frame and formula).

You will find the changes in the last version of the svn (1.7.2). Let me
know if it works for you or if you have additional questions.

Best regards,
Alejandro Reyes


> Hi,
>    Would someone be kind enough to clarify the two questions raised about DEXSeq in my earlier email below?
>
>    Would attaching any further information help? I am working with DEXSeq 1.6.0, and here are the model frame and model matrices if it makes it easier. Thanks!
>
>> mf
>          sample exon sizeFactor condition   batch dispersion count
> 1 Untr_biorep1 E001  1.0517803      Untr biorep1 0.13523549    25
> 2 Untr_biorep1 E002  1.0517803      Untr biorep1 0.01656906   708
> 3  LPS_biorep1 E001  1.0010474       LPS biorep1 0.13523549    20
> 4  LPS_biorep1 E002  1.0010474       LPS biorep1 0.01656906   442
> 5 Untr_biorep2 E001  0.8920483      Untr biorep2 0.13523549    26
> 6 Untr_biorep2 E002  0.8920483      Untr biorep2 0.01656906   947
> 7  LPS_biorep2 E001  1.1090401       LPS biorep2 0.13523549    25
> 8  LPS_biorep2 E002  1.1090401       LPS biorep2 0.01656906   884
>
>> mm
>    (Intercept) sampleLPS_biorep2 sampleUntr_biorep1 sampleUntr_biorep2
> 1           1                 0                  1                  0
> 2           1                 0                  1                  0
> 3           1                 0                  0                  0
> 4           1                 0                  0                  0
> 5           1                 0                  0                  1
> 6           1                 0                  0                  1
> 7           1                 1                  0                  0
> 8           1                 1                  0                  0
>    conditionUntr batchbiorep2 exonE002 conditionUntr:exonE002
> 1             1            0        0                      0
> 2             1            0        1                      1
> 3             0            0        0                      0
> 4             0            0        1                      0
> 5             1            1        0                      0
> 6             1            1        1                      1
> 7             0            1        0                      0
> 8             0            1        1                      0
>    batchbiorep2:exonE002
> 1                     0
> 2                     0
> 3                     0
> 4                     0
> 5                     0
> 6                     1
> 7                     0
> 8                     1
> attr(,"assign")
> [1] 0 1 1 1 2 3 4 5 6
> attr(,"contrasts")
> attr(,"contrasts")$sample
> [1] "contr.treatment"
>
> attr(,"contrasts")$condition
> [1] "contr.treatment"
>
> attr(,"contrasts")$batch
> [1] "contr.treatment"
>
> attr(,"contrasts")$exon
> [1] "contr.treatment"
>
>
>
>
>
>
>
> From: Narayanan, Manikandan (NIH/NIAID) [E]
> Sent: Thursday, May 16, 2013 11:14 AM
> To: bioconductor at r-project.org
> Subject: DEXSeq on two-exon genes: how to specify a formula without redundant terms
>
> Hi DEXSeq users/developers,
>    I have used DEXSeq successfuly for genes with many exons and really like the diagnostic/visualization plots that come with it. Recently though, for genes with two testable exons, I am getting the "Underdetermined model; cannot estimate dispersions." error.
>
>    I figure this is due to redundant terms in my formula as shown in PS below. So my questions are:
>
> 1) Is there a way to specify the formula  count ~ sample + (condition + batch) * exon so that redundant terms 'condition + batch' are removed?
>
> 2) If not, is it safe to change ncol(mm) to qr(mm)$rank (i.e., rank of model matrix to remove redundant terms) in this piece of code in estimateExonDispersionsForModelFrame:
>      if (nrow(mm) <= ncol(mm))
>          stop("Underdetermined model; cannot estimate dispersions. Maybe replicates have not been properly specified.")
>
> Would changing the code this way violate any assumptions of the DEXSeq model?
>
>
> Thank you,
> Mani
>
>
> PS: # condition + batch terms are redundant as sample term is already present!
>> formulaDispersion
> count ~ sample + (condition + batch) * exon
>
>> design(ecs)
>                   condition       batch
> Untr_biorep1      Untr     biorep1
> LPS_biorep1        LPS     biorep1
> Untr_biorep2      Untr     biorep2
> LPS_biorep2        LPS     biorep2
>
>> colnames(model.matrix(formulaDispersion, mf))
> [1] "(Intercept)"            "sampleLPS_biorep2"      "sampleUntr_biorep1"
> [4] "sampleUntr_biorep2"     "conditionUntr"          "batchbiorep2"
> [7] "exonE002"               "conditionUntr:exonE002" "batchbiorep2:exonE002"
>
>       [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor




More information about the Bioconductor mailing list