[Bioc-devel] -which() pitfall in subsetting exressions
Wolfgang Huber
whuber at embl.de
Sun May 26 15:59:45 CEST 2013
[bioc-devel added and subject header changed]
Dear List & Developers
Mani's post reminded me of the pitfall posed by using expressions such as
x = x[ -which(someExpression) ]
which tends to have unintended consequences in particular if all elements of 'someExpression' are FALSE. Compare:
x=1:10
> x[ -which(x==17) ]
integer(0)
> x[ x!=17 ]
[1] 1 2 3 4 5 6 7 8 9 10
A quick search in Bioconductor source code found this idiom in multiple packages, including:
IRanges Biostrings beadarray BeadDataPackR BioNet affycoretools annotate copynumber CRImage DiffBind flowStats geNetClassifier GGtools maigesPack maskBAD MineICA MmPalateMiRNA Mulcom nem phyloseq plgem RBGL girafe Ringo RIPSeeker rnaSeqMap SomatiCA Starr VanillaICE
I cannot tell whether this poses an actual bug threat in any of these cases, but developers might want to double-check.
Best wishes
Wolfgang
On May 24, 2013, at 5:08 pm, "Narayanan, Manikandan (NIH/NIAID) [E]" <manikandan.narayanan at nih.gov> wrote:
> 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
>
>
> _______________________________________________
> 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 Bioc-devel
mailing list