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

Martin Morgan martin.morgan at roswellpark.org
Thu Mar 17 15:35:03 CET 2016



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



More information about the Bioc-devel mailing list