[Bioc-devel] Making GenomicAlignments::readGAlignmentPairs() fail fast if given bad seqlevels in `which`

[Martin Morgan]

On 03/17/2016 10:35 AM, Martin Morgan wrote:
>
>
> On 03/17/2016 10:29 AM, Peter Hickey wrote:
>> Thanks, Aaron. I implemented a similar workaround, but I think it
>> would be nice to have in the core Bioconductor implementation. I had a
>> quick poke around GenomicAlignments::readGAlignmentPairs(), however,
>> but it looked like I'd have to learn a bit too much about the
>> underlying Rsamtools::scanBam() in order to implement a quick fix.
>>
>
> The error does come from scanBam
>
>  > bf = BamFile(system.file(package="Rsamtools", "extdata", "ex1.bam"))
>  > gr <- GRanges(paste0("seq", 3), IRanges(1, 1000))
>  > param <- ScanBamParam(which=gr, what="qname")
>  > scanBam(bf, param=param)
> Error in value[[3L]](cond) : 'scanBam' failed:
> record: 0
> error: 0
> file:
> /home/mtmorgan/R/x86_64-pc-linux-gnu-library/3.3/Rsamtools/extdata/ex1.bam
> index:
> /home/mtmorgan/R/x86_64-pc-linux-gnu-library/3.3/Rsamtools/extdata/ex1.bam.bai
>
> In addition: Warning message:
> In doTryCatch(return(expr), name, parentenv, handler) :
> space 'seq3' not in BAM header
>
> and it seems like this is an easy fix. Look for it later today.

in svn at 1.23.5, and in Bioc-devel hopefully after tonight's build.

Martin

>
> @Aaron -- better to use seqlevels(BamFile(.)) rather than parsing the
> target from scanBamHeader), both to avoid code duplication and to
> benefit from whatever bugs (e.g., when bam files have duplicate header
> entries) are reported.
>
>  > selectMethod(seqinfo, "BamFile")
> Method Definition:
>
> function (x)
> {
>      h <- scanBamHeader(x, what = "targets")[["targets"]]
>      h <- h[!(duplicated(h) & duplicated(names(h)))]
>      Seqinfo(names(h), unname(h))
> }
> <environment: namespace:Rsamtools>
>
> Signatures:
>          x
> target  "BamFile"
> defined "BamFile"
>
> Martin
>
>>>
>>> Hi Peter,
>>>
>>> I had the same problem a while ago and solved it by first reading
>>> only the
>>> header of the BAM file, extracting the chromosomes that are available
>>> and
>>> generating a warning for all given chromosomes that are not
>>> available. That
>>> worked for my purposes. I have implemented this in a function (
>>> https://github.com/ataudt/aneufinder/blob/master/R/importReads.R)
>>>
>>> Aaron
>>>
>>>          [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> Bioc-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/bioc-devel
>>
>
>
> This email message may contain legally privileged and/or confidential
> information.  If you are not the intended recipient(s), or the employee
> or agent responsible for the delivery of this message to the intended
> recipient(s), you are hereby notified that any disclosure, copying,
> distribution, or use of this email message is prohibited.  If you have
> received this message in error, please notify the sender immediately by
> e-mail and delete this email message from your computer. Thank you.
>
> _______________________________________________
> Bioc-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/bioc-devel


This email message may contain legally privileged and/or confidential information.  If you are not the intended recipient(s), or the employee or agent responsible for the delivery of this message to the intended recipient(s), you are hereby notified that any disclosure, copying, distribution, or use of this email message is prohibited.  If you have received this message in error, please notify the sender immediately by e-mail and delete this email message from your computer. Thank you.