<p>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 <code>scanBam</code> in different windows (<code>which</code> parameter): the first call to <code>scanBam</code> succeeds but subsequent calls return zero-length results.</p>
<p>By contrast, if I don’t specify a coordinate window and instead use <code>yieldSize</code> when opening the BAM file, I can successfully read records in a loop:</p>
<pre><code>library(Rsamtools)
bamfile = system.file('extdata', 'ex1.bam', package = 'Rsamtools')
process = function (reads) message(length(reads$seq))
</code></pre>

<pre><code>bam = open(BamFile(bamfile, yieldSize = 1000L))
param = ScanBamParam(what = 'seq')

repeat {
    reads = scanBam(bam, param = param)[[1L]]
    if (length(reads$seq) == 0L) break
    process(reads)
}
</code></pre>
<p>… where <code>process</code> performs the read processing (unimportant for the example).</p>
<p>Now I’m trying to use the following code to work in genomic coordinate windows instead, and it fails as explained above:</p>
<pre><code>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]]
    process(reads)
})
</code></pre>
<p>This only returns data from the first iteration, not any subsequent iterations.</p>
<p>Conversely, it works if I close and re-open the BAM file in each iteration (i.e. before each <code>scanBam</code> call) — but that’s undesirable: it carries a <strong>substantial</strong> overhead (running on an <em>empty</em> BAM file with 16 chromosomes in its header takes several <em>hours</em>).</p>
<p>Since the <code>process</code> function is nontrivial and written in R, I’m stuck with R for now. And there’s no reason why this <em>shouldn’t</em> work in R.</p>
<p>Am I overlooking something? Is this a bug in {Rsamtools}? Is it “by design”?</p>

CC BY-SA 4.0
2021-07-08T08:57:39.643Z
@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](https://meta.stackexchange.com/q/19190/1968), 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).
Konrad Rudolph
CC BY-SA 4.0
2021-07-08T10:53:25.270Z
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.
gringer
CC BY-SA 4.0
2021-10-20T07:49:19.177Z
Are you working under linux? and do you absolutely have to do it in R?
user324810
CC BY-SA 4.0
2021-10-20T08:12:06.030Z
@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.
Konrad Rudolph
<r><sam>
How can I read multiple different regions from a BAM file in R?
# How can I read multiple different regions from a BAM file in R?

(Asked by: username_0 on 2018-10-04)

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:
library(Rsamtools)
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
    process(reads)
}

… 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]]
    process(reads)
})

This only returns data from the first iteration, not any subsequent iterations.
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).
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”?

**Tags:** r, sam
> **Comments:**
> **username_0**: @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](https://meta.stackexchange.com/q/19190/1968), 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).  
> **username_1**: 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.  
> **username_2**: Are you working under linux? and do you absolutely have to do it in R?  
> **username_0**: @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.  

---

## 0 answers

