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

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.

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]]

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.

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

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.