[Bioc-devel] Slow performance on scanBam

I've been dealing with some small bam files with millions of reference
sequences leading to monster headers. As one might guess, this leads to
pretty bad performance when calling scanBam.

Right now, it takes a bit (27MB bam file, 16k alignments, 2.5 million
reference sequences in the reference fasta file):

scanBam("sim.fasta-L27-ma8-mp6-rfg5_2-rdg3_1.bam")

    user  system elapsed
186.264   0.528 186.934

I've traced it down to scanBamHeader and seqinfo-BamFile, the problematic
code is in scanBamHeader which processes the entire header when all seqinfo
needs is the `targets` portion of the list. Additionally, the
order(rankSeqLevels(.)) doesn't scale either. So I've replaced this as
well. I've changed the body of seqinfo-BamFile from:

     h <- scanBamHeader(x)[["targets"]]
     o <- order(rankSeqlevels(names(h)))
     Seqinfo(names(h)[o], unname(h)[o])

to this:

     if (!isOpen(x)) {
         open(x)
         on.exit(close(x))
     }
     h <- .Call(.read_bamfile_header, .extptr(x))$targets
     Seqinfo(names(h), unname(h))

We then get:

scanBam("sim.fasta-L27-ma8-mp6-rfg5_2-rdg3_1.bam")

    user  system elapsed
  14.780   0.360  15.158

Which is still pretty slow for how small these files are, but I can
probably live with that.

Two questions:
-- do we need the: order(rankSeqlevels(names(h))) bit? that does change the
return value, but I can certainly live with the ordering in the file.
-- what else am I missing?

I can send a patch if need be, but would like to hear from the cognoscenti
first if there is a "built-in" way to avoid this overhead.

thanks, jim

	[[alternative HTML version deleted]]

Hi Martin, thanks for the quick response. The data is certainly shareable. Here
is a link to a bam + bai + sam file that I have been using for benchmarking:
https://www.dropbox.com/s/eat31mnmmco1zoh/example-bam.tar.bz2

There is a method in SAM to elide the reference names from a header, but I think
they are just shunted to another file so I gave up on that track. Since I only
end up aligning to a small fraction of the transcripts, I might be able to
post-process the file, but it would be best to process as-is.

thanks again, jim



On Tue, May 13, 2014 at 5:16 AM, Martin Morgan <mtmorgan at fhcrc.org
<mailto:mtmorgan at fhcrc.org>> wrote:

    Hi James -- I don't think there's anything in existence to make this easier,
    but I'll expose something in the next 24 hours; is your data shareable?
    There might be deeper things to be done for processing this
    small-but-numerous style data.

    Martin


    On 05/12/2014 05:32 PM, James Bullard wrote:

        I've been dealing with some small bam files with millions of reference
        sequences leading to monster headers. As one might guess, this leads to
        pretty bad performance when calling scanBam.

        Right now, it takes a bit (27MB bam file, 16k alignments, 2.5 million
        reference sequences in the reference fasta file):

            scanBam("sim.fasta-L27-ma8-__mp6-rfg5_2-rdg3_1.bam")

             user  system elapsed
        186.264   0.528 186.934

        I've traced it down to scanBamHeader and seqinfo-BamFile, the problematic
        code is in scanBamHeader which processes the entire header when all seqinfo
        needs is the `targets` portion of the list. Additionally, the
        order(rankSeqLevels(.)) doesn't scale either. So I've replaced this as
        well. I've changed the body of seqinfo-BamFile from:

              h <- scanBamHeader(x)[["targets"]]
              o <- order(rankSeqlevels(names(h)))
              Seqinfo(names(h)[o], unname(h)[o])

        to this:

              if (!isOpen(x)) {
                  open(x)
                  on.exit(close(x))
              }
              h <- .Call(.read_bamfile_header, .extptr(x))$targets
              Seqinfo(names(h), unname(h))

        We then get:

            scanBam("sim.fasta-L27-ma8-__mp6-rfg5_2-rdg3_1.bam")

             user  system elapsed
           14.780   0.360  15.158

        Which is still pretty slow for how small these files are, but I can
        probably live with that.

        Two questions:
        -- do we need the: order(rankSeqlevels(names(h))) bit? that does change the
        return value, but I can certainly live with the ordering in the file.
        -- what else am I missing?

        I can send a patch if need be, but would like to hear from the cognoscenti
        first if there is a "built-in" way to avoid this overhead.

        thanks, jim

                 [[alternative HTML version deleted]]