6
$\begingroup$

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)
})

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).

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

$\endgroup$
  • $\begingroup$ Have you tried using the plyranges? They have a specific function for rolling a fixed window, can be used on a bam file. $\endgroup$ – llrs Oct 9 '18 at 14:43
  • $\begingroup$ What is the end goal? To get all of the reads aligning to a specific window? Also, are you open to using solutions that are not in R? (python, unix tools, et cetera) $\endgroup$ – conchoecia Oct 9 '18 at 15:53
  • $\begingroup$ @conchoecia 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. $\endgroup$ – Konrad Rudolph Oct 10 '18 at 8:25
  • $\begingroup$ @Llopis Correct me if I’m wrong but doesn’t plyranges need the data loaded into memory as a whole? Unfortunately the memory overhead of that may be prohibitive. I’d like to retain the streaming access. That said, I will try that. $\endgroup$ – Konrad Rudolph Oct 10 '18 at 8:27
  • 1
    $\begingroup$ I never used it, I'm not into the genomic world.The vignette doesn't say anything, but in the manual it says "Reading a BAM file is deferred until an action" is done. So it might only load what is needed to perform the tasks $\endgroup$ – llrs Oct 10 '18 at 9:49

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.