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]]
_______________________________________________
[email protected] 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.
_______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/bioc-devel
