I’m trying to process BAM records (of a coordinate-sorted, indexed BAM file) successively in fixed-size, non-overlapping, sliding genome coordinate windows. Unfortunately I am unable to read records via scanBam in different windows (which parameter): the first call to scanBam succeeds but subsequent calls return zero-length results.

By contrast, if I don’t specify a coordinate window and instead use yieldSize when opening the BAM file, I can successfully read records in a loop:

bamfile = system.file('extdata', 'ex1.bam', package = 'Rsamtools')
process = function (reads) message(length(reads$seq))
bam = open(BamFile(bamfile, yieldSize = 1000L))
param = ScanBamParam(what = 'seq')

repeat {
    reads = scanBam(bam, param = param)[[1L]]
    if (length(reads$seq) == 0L) break

… where process performs the read processing (unimportant for the example).

Now I’m trying to use the following code to work in genomic coordinate windows instead, and it fails as explained above:

windowsize = 1000L

bam = open(BamFile(bamfile))
chr = seqlevels(bam)[1L] # Only do one chromosome for now
chrlengths = seqlengths(bam)

result = lapply(seq(1L, chrlengths[chr] - windowsize, by = windowsize), function (pos) {
    which = setNames(RangesList(IRanges(pos, pos + windowsize)), chr)
    param = ScanBamParam(what = 'seq', which = which)
    reads = scanBam(bam, param = param)[[1L]]

Conversely, it works if I close and re-open the BAM file in each iteration (i.e. before each scanBam call) — but that’s undesirable: it carries a substantial overhead (running on an empty BAM file with 16 chromosomes in its header takes several hours).

It's only returning data from the first iteration, not any subsequent iterations, unless I close the BAM file between scanBam calls.

Since the process function is nontrivial and written in R, I’m stuck with R for now. And there’s no reason why this shouldn’t work in R.

Am I overlooking something? Is this a bug in {Rsamtools}? Is it “by design”?

  • $\begingroup$ @gringer I’ve noticed that you’ve added the tag into the question title — why? On Stack Overflow there’s a broad consensus against redundant mentions of the tag in the question title, and in fact I frequently (several times a day, on average) edit to remove such mentions from titles on Stack Overflow (most useless taggings are more egregious, but removing “… in X” is a fair share of the edits I make). $\endgroup$ Jul 8 at 8:57
  • 1
    $\begingroup$ Following the same logic, why not remove BAM from the title? I don't think that "in R" is redundant here. In this case I considered that there were multiple valid Bioinformatics questions that could begin with "How can I read multiple different regions from a BAM file", even with the "R" tag. In this case it seems like you specifically want an R solution, which is why I added it in. $\endgroup$
    – gringer
    Jul 8 at 10:53
  • $\begingroup$ Are you working under linux? and do you absolutely have to do it in R? $\endgroup$
    – user324810
    Oct 20 at 7:49
  • $\begingroup$ @user324810 I mean, this question is three years old. And at any rate, regardless of whether somebody “absolutely” has to do this in R, it should work in R. But for posterity: yes, I absolutely had to do this in R since the process placeholder in my case was a piece of nontrivial logic developed in R over the span of months (in fact, the final implementation was in C++ anyway, the R code was my research prototype for faster iteration). I ended up processing fixed-sized chunks rather than scanning over genomic coordinate windows. $\endgroup$ Oct 20 at 8:12
  • $\begingroup$ Didn't notice the time. I filtered the questions by unanswered and I saw it hanging there. Sorry for re-opening it. In any case, glad you could a solution to your problem. I would have proposed to simply execute a bash command with samtools view if you are under Linux OS. $\endgroup$
    – user324810
    Oct 20 at 8:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.