[Bioc-devel] GenomicFiles: chunking

Valerie Obenchain vobencha at fhcrc.org
Mon Oct 27 19:35:46 CET 2014

Hi Kasper and Mike,

I've added 2 new functions to GenomicFiles and deprecated the old 
classes. The vignette now has a graphic which (hopefully) clarifies the 
MAP / REDUCE mechanics of the different functions.

Below is some performance testing for the new functions and answers to 
leftover questions from previous emails.

Major changes to GenomicFiles (in devel):

- *FileViews classes have been deprecated:

The idea is to use the GenomicFiles class to hold any type of file be it 
BAM, BigWig, character vector etc. instead of having name-specific 
classes like BigWigFileViews. Currently GenomicFiles does not inherit 
from SummarizedExperiment but it may in the future.

- Add reduceFiles() and reduceRanges():

These functions pass all ranges or files to MAP vs the lapply approach 
taken in reduceByFiles() and reduceByRanges().

(1) Performance:

When testing with reduceByFile() you noted "GenomicFiles" is 10-20x 
slower than the straightforward approach". You also noted this was 
probably because of the lapply over all ranges - true. (Most likely 
there was overhead in creating the SE object as well.) With the new 
reduceFiles(), passing all ranges at once, we see performance very 
similar to that of the 'by hand' approach.

In the test code I've used Bam instead of BigWig. Both test functions 
output lists, have comparable MAP and REDUCE steps etc.

I used 5 files ('fls') and a granges ('grs') of length 100.
 > length(grs)
[1] 100

 > sum(width(grs))
[1] 1000000

FUN1 is the 'by hand' version. These results are similar to what you 
saw, not quite a 4x difference between 10 and 100 ranges.

 >> microbenchmark(FUN1(grs[1:10], fls), FUN1(grs[1:100], fls), times=10)
 > Unit: seconds
 >                   expr      min       lq     mean   median       uq 
 >   FUN1(grs[1:10], fls) 1.177858 1.190239 1.206127 1.201331 1.222298 
 >  FUN1(grs[1:100], fls) 4.145503 4.163404 4.249619 4.208486 4.278463 
 >  neval
 >     10
 >     10

FUN2 is the reduceFiles() approach and the results are very similar to FUN1.

 >> microbenchmark(FUN2(grs[1:10], fls), FUN2(grs[1:100], fls), times=10)
 > Unit: seconds
 >                   expr      min       lq     mean   median       uq 
 >   FUN2(grs[1:10], fls) 1.242767 1.251188 1.257531 1.253154 1.267655 
 >  FUN2(grs[1:100], fls) 4.251010 4.340061 4.390290 4.361007 4.384064 
 >  neval
 >     10
 >     10

(2) Request for "chunking of the mapping of ranges":

For now we decided not to add a 'yieldSize' argument for chunking. There 
are 2 approaches to chunking through ranges *within* the same file. In 
both cases the user spits the ranges, either before calling the function 
or in the MAP step.

i) reduceByFile() with a GRangesList:

The user provides a GRangesList as the 'ranges' arg. On each worker, 
lapply applies MAP to the one files and all elements of the GRangesList.

ii) reduceFiles() with a MAP that handles chunking:

The user split ranges in MAP and uses lapply or another loop to iterate. 
For example,

MAP <- function(range, file, ...) {
     lst = split(range, someFactor)
     someFUN = function(rngs, file, ...) do something
     lapply(lst, FUN=someFun, ...)

The same ideas apply for chunking though ranges *across* files with 
reduceByRange() and reduceRanges().

iii) reduceByRange() with a GRangesList:

Mike has a good example here:

iv) reduceRanges():

'ranges' should be a GRangesList. The MAP step will operate on an 
element of the GRangesList and all files. Unless you want to operate on 
all files at once I'd use reduceByRange() instead.

(3) Return objects have different shape:

Previous question:

Why does the return object of reduceByFile vs reduceByRange (with
summarize = FALSE) different?  I understand why internally you have
different nesting schemes (files and ranges) for the two functions, but 
it is not clear to me that it is desirable to have the return object 
depend on how the computation was done.

reduceByFile() and reduceFiles() output a list the same length as the 
number of files while reduceByRange() and reduceRanges() output a list 
the same length as the number of ranges.

Reduction is different depending on which function is chosen; data are 
collapsed either within a file or across files. When REDUCE does 
something substantial the output are not equivalent.

While it's possible to get the same result (REDUCE simply unlists or 
isn't used), the two approaches were not intended to be equivalent ways 
of arriving at the same end. The idea was that the user had a specific 
use case in mind - they either wanted to collapse the data across or 
within files.

(4) return object from coverage(BigWigFileViews):

Previous comment:

coverage(BigWigFileViews) returns a "wrong" assay object in my opinion,
Specifically, each (i,j) entry in the object is an RleList with a single
element with a name equal to the seqnames of the i'th entry in the query
GRanges.  To me, this extra nestedness is unnecessary; I would have
expected an Rle instead of an RleList with 1 element.

The return value from coverage(x) is an RleList with one coverage vector 
per seqlevel in 'x'. Even if there is only one seqlevel, the result 
still comes back as an RleList. This is just the default behavior.

(5) separate the 'read' function from the MAP step

Previous comment:

Also, something completely different, it seems like it would be 
convenient for stuff like BigWigFileViews to not have to actually parse 
the file in the MAP step.  Somehow I would envision some kind of reading 
function, stored inside the object, which just returns an Rle when I ask 
for a (range, file).  Perhaps this is better left for later.

The current approach for the reduce* functions is for MAP to both 
extract and manipulate data. The idea of separating the extraction step 
is actually implemented in reduceByYield(). (This function used to be 
yieldReduce() in Rsamtools in past releases.) For reduceByYield() teh 
user must specify YIELD (a reader function), MAP, REDUCE and DONE 
(criteria to stop iteration).

I'm not sure what is best here. I thought the many-argument approach of 
reduceByYield() was possibly confusing or burdensome and so didn't use 
it in the other GenomicFiles functions. Maybe it's not confusing but 
instead makes the individual steps more clear. What do you think,

- Should the reader function be separate from the MAP? What are the 

- Should READER, MAP, REDUCE be stored inside the GenomicFiles object or 
supplied as arguments to the functions?

(6) unnamed assay in SummarizedExperiment

Previous comment:

The return object of reduceByRange / reduceByFile with summarize = TRUE 
is a SummarizedExperiment with an unnamed assay.  I was surprised to see 
that this is even possible.

There is no default name for SummarizedExperiment in general. I've named 
the assay 'data' for lack of a better term. We could also go with 
'reducedData' or another suggestion.

Thanks for the feedback.


On 10/01/2014 08:30 AM, Michael Love wrote:
> hi Kasper,
> For a concrete example, I posted a R and Rout file here:
> https://gist.github.com/mikelove/deaff999984dc75f125d
> Things to note: 'ranges' is a GRangesList, I cbind() the numeric
> vectors in the REDUCE, and then rbind() the final list to get the
> desired matrix.
> Other than the weird column name 'init', does this give you what you want?
> best,
> Mike
> On Tue, Sep 30, 2014 at 2:08 PM, Michael Love
> <michaelisaiahlove at gmail.com> wrote:
>> hi Kasper and Valerie,
>> In Kasper's original email:
>> "I would like to be able to write a MAP function which takes
>>    ranges, file
>> instead of just
>>    range, file
>> And then chunk over say 1,000s of ranges. I could then have an
>> argument to reduceByXX called something like rangeSize, which is kind
>> of yieldSize."
>> I was thinking, for your BigWig example, we get part of the way just
>> by splitting the ranges into a GRangesList of your desired chunk size.
>> Then the parallel iteration is over chunks of GRanges. Within your MAP
>> function:
>> import(file, which = ranges, as = "Rle", format = "bw")[ranges]
>> returns an RleList, and calling mean() on this gives a numeric vector
>> of the mean of the coverage for each range.
>> Then the only work needed is how to package the result into something
>> reasonable. What you would get now is a list (length of GRangesList)
>> of lists (length of files) of vectors (length of GRanges).
>> best,
>> Mike
>> On Mon, Sep 29, 2014 at 7:09 PM, Valerie Obenchain <vobencha at fhcrc.org> wrote:
>>> Hi Kasper,
>>> The reduceBy* functions were intended to combine data across ranges or
>>> files. If we have 10 ranges and 3 files you can think of it as a 10 x 3 grid
>>> where we'll have 30 queries and therefore 30 pieces of information that can
>>> be combined across different dimensions.
>>> The request to 'processes ranges in one operation' is a good suggestion and,
>>> as you say, may even be the more common use case.
>>> Some action items and explanations:
>>> 1)  treat ranges / files as a group
>>> I'll add functions to treat all ranges / files as a group; essentially no
>>> REDUCER other than concatenation. Chunking (optional) will occur as defined
>>> by 'yieldSize' in the BamFile.
>>> 2)  class clean-up
>>> The GenomicFileViews class was the first iteration. It was overly
>>> complicated and the complication didn't gain us much. In my mind the
>>> GenomicFiles class is a striped down version should replace the *FileViews
>>> classes. I plan to deprecate the *FileViews classes.
>>> Ideally the class(s) in GenomicFiles would inherit from
>>> SummarizedExperiment. A stand-alone package for SummarizedExperiment is in
>>> the works for the near future. Until those plans are finalized I'd rather
>>> not change the inheritance structure in GenomicFiles.
>>> 3) pack / unpack
>>> I experimented with this an came to the conclusion that packing / unpacking
>>> should be left to the user vs being done behind the scenes with reduceBy*.
>>> The primary difficulty is that you don't have pre-knowledge of the output
>>> format of the user's MAPPER. If MAPPER outputs an Rle, unpacking may be
>>> straightforward. If MAPPER outputs a NumericList or matrix or vector with no
>>> genomic coordinates then things are more complicated.
>>> I'm open if others have suggestions / prototypes.
>>> 4) reduceByYield
>>> This was intended for chunking through ranges in a single file. You can
>>> imagine using bplapply() over files where each file is chunked through with
>>> reduceByYield(). All ranges are chunked by 'yieldSize' defined in the
>>> BamFile unless a 'param' dictates a subset of ranges.
>>> 5) additional tidy
>>> I'll add a name to the 'assays' when summarize=TRUE, make sure return types
>>> are consistent when summarize=FALSE, etc.
>>> Thanks for the testing and feedback.
>>> Valerie
>>> On 09/29/2014 07:18 AM, Michael Love wrote:
>>>> On Mon, Sep 29, 2014 at 9:09 AM, Kasper Daniel Hansen
>>>> <kasperdanielhansen at gmail.com> wrote:
>>>>> I don't fully understand "the use case for reducing by range is when the
>>>>> entire dataset won't fit into memory".  The basic assumption of these
>>>>> functions (as far as I can see) is that the output data fits in memory.
>>>>> What may not fit in memory is various earlier "iterations" of the data.
>>>>> For
>>>>> example, in my use case, if I just read in all the data in all the ranges
>>>>> in
>>>>> all the samples it is basically Rle's across 450MB times 38 files, which
>>>>> is
>>>>> not small.  What fits in memory is smaller chunks of this; that is true
>>>>> for
>>>>> every application.
>>>> I was unclear. I meant that, in approach1, you have an object,
>>>> all.Rle, which contains Rles for every range over every file. Can you
>>>> actually run this approach on the full dataset?
>>>>> Reducing by range (or file) only makes sense when the final output
>>>>> includes
>>>>> one entity for several ranges/files ... right?  So I don't see how reduce
>>>>> would help me.
>>>> Yes, I think we agree. This is not a good use case for reduce by range
>>>> as now implemented.
>>>> This is a use case which would benefit from the user-facing function
>>>> internally calling pack()/unpack() to reduce the number of import()
>>>> calls, and then in the end giving back the mean coverage over the
>>>> input ranges. I want this too.
>>>> https://github.com/Bioconductor/GenomicFileViews/issues/2#issuecomment-32625456
>>>> (link to the old github repo, the new github repo is named GenomicFiles)
>>>>> As I see the pack()/unpack() paradigm, it just re-orders the query ranges
>>>>> (which is super nice and matters a lot for speed for some applications).
>>>>> As
>>>>> I understand the code (and my understanding is developing) we need an
>>>>> extra
>>>>> layer to support processing multiple ranges in one operation.
>>>>> I am happy to help apart from complaining.
>>>>> Best,
>>>>> Kasper
>>>>> On Mon, Sep 29, 2014 at 8:55 AM, Michael Love
>>>>> <michaelisaiahlove at gmail.com>
>>>>> wrote:
>>>>>> Thanks for checking it out and benchmarking. We should be more clear
>>>>>> in the docs that the use case for reducing by range is when the entire
>>>>>> dataset won't fit into memory. Also, we had some discussion and
>>>>>> Valerie had written up methods for packing up the ranges supplied by
>>>>>> the user into a better form for querying files. In your case it would
>>>>>> have packed many ranges together, to reduce the number of import()
>>>>>> calls like your naive approach. See the pack/unpack functions, which
>>>>>> are not in the vignette but are in the man pages. If I remember,
>>>>>> writing code to unpack() the result was not so simple, and development
>>>>>> of these functions was set aside for the moment.
>>>>>> Mike
>>>>>> On Sun, Sep 28, 2014 at 10:49 PM, Kasper Daniel Hansen
>>>>>> <kasperdanielhansen at gmail.com> wrote:
>>>>>>> I am testing GenomicFiles.
>>>>>>> My use case: I have 231k ranges of average width 1.9kb and total width
>>>>>>> 442
>>>>>>> MB.  I also have 38 BigWig files.  I want to compute the average
>>>>>>> coverage
>>>>>>> of the 38 BigWig files inside each range.  This is similar to wanting
>>>>>>> to
>>>>>>> get coverage of - say - all promoters in the human genome.
>>>>>>> My data is residing on a file server which is connected to the compute
>>>>>>> node
>>>>>>> through ethernet, so basically I have very slow file access.  Also. the
>>>>>>> BigWig files are (perhaps unusually) big: ~2GB / piece.
>>>>>>> Below I have two approaches, one using basically straightforward code
>>>>>>> and
>>>>>>> the other using GenomicFiles (examples are 10 / 100 ranges on 5 files).
>>>>>>> Basically GenomicFiles is 10-20x slower than the straightforward
>>>>>>> approach.
>>>>>>> This is most likely because reduceByRange/reduceByFile processes each
>>>>>>> (range, file) separately.
>>>>>>> It seems naturally (to me) to allow some chunking of the mapping of the
>>>>>>> ranges.  My naive approach is fast (I assume) because I read multiple
>>>>>>> ranges through one import() command.  I know I have random access to
>>>>>>> the
>>>>>>> BigWig files, but there must still be some overhead of seeking and
>>>>>>> perhaps
>>>>>>> more importantly instantiating/moving stuff back and forth to R.  So
>>>>>>> basically I would like to be able to write a MAP function which takes
>>>>>>>     ranges, file
>>>>>>> instead of just
>>>>>>>     range, file
>>>>>>> And then chunk over say 1,000s of ranges.  I could then have an
>>>>>>> argument
>>>>>>> to
>>>>>>> reduceByXX called something like rangeSize, which is kind of yieldSize.
>>>>>>> Perhaps this is what is intended for the reduceByYield on BigWig files?
>>>>>>> In a way, this is what is done in the vignette with the
>>>>>>> coverage(BAMFILE)
>>>>>>> example where tileGenome is essentially constructed by the user to
>>>>>>> chunk
>>>>>>> the coverage computation.
>>>>>>> I think the example of a set of regions I am querying on the files,
>>>>>>> will
>>>>>>> be
>>>>>>> an extremely common usecase going forward. The raw data for all the
>>>>>>> regions
>>>>>>> together is "too big" but I do a computation on each region to reduce
>>>>>>> the
>>>>>>> size.  In this situation, all the focus is on the MAP step.  I see the
>>>>>>> need
>>>>>>> for REDUCE in the case of the t-test example in the vignette, where the
>>>>>>> return object is a single "thing" for each base.  But in general, I
>>>>>>> think
>>>>>>> we will use these file structures a lot to construct things without
>>>>>>> REDUCE
>>>>>>> (across neither files nor ranges).
>>>>>>> Also, something completely different, it seems like it would be
>>>>>>> convenient
>>>>>>> for stuff like BigWigFileViews to not have to actually parse the file
>>>>>>> in
>>>>>>> the MAP step.  Somehow I would envision some kind of reading function,
>>>>>>> stored inside the object, which just returns an Rle when I ask for a
>>>>>>> (range, file).  Perhaps this is better left for later.
>>>>>>> Best,
>>>>>>> Kasper
>>>>>>> Examples
>>>>>>> approach1 <- function(ranges, files) {
>>>>>>>       ## By hand
>>>>>>>       all.Rle <- lapply(files, function(file) {
>>>>>>>           rle <- import(file, as = "Rle", which = ranges, format =
>>>>>>> "bw")[ranges]
>>>>>>>           rle
>>>>>>>       })
>>>>>>>       print(object.size(all.Rle), units = "auto")
>>>>>>>       mat <- do.call(cbind, lapply(all.Rle, function(xx) {
>>>>>>>           sapply(xx, mean)
>>>>>>>       }))
>>>>>>>       invisible(mat)
>>>>>>> }
>>>>>>> system.time({
>>>>>>> mat1 <- approach1(all.grs[1:10], all.files[1:5])
>>>>>>> })
>>>>>>> 160.9 Kb
>>>>>>>       user  system elapsed
>>>>>>>      1.109   0.001   1.113
>>>>>>> system.time({
>>>>>>> mat1 <- approach1(all.grs[1:100], all.files[1:5])
>>>>>>> }) # less than 4x slower than previous call
>>>>>>> 3 Mb
>>>>>>>       user  system elapsed
>>>>>>>      4.101   0.019   4.291
>>>>>>> approach2 <- function(ranges, files) {
>>>>>>>       gf <- GenomicFiles(rowData = ranges, files = files)
>>>>>>>       MAPPER <- function(range, file, ....) {
>>>>>>>           library(rtracklayer)
>>>>>>>           rle <- import(file, which = range, as = "Rle", format =
>>>>>>> "bw")[range]
>>>>>>>           mean(rle)
>>>>>>>       }
>>>>>>>       sExp <- reduceByRange(gf, MAP = MAPPER, summarize = TRUE, BPPARAM
>>>>>>> =
>>>>>>> SerialParam())
>>>>>>>       sExp
>>>>>>> }
>>>>>>> system.time({
>>>>>>> mat2 <- approach2(all.grs[1:10], all.files[1:5])
>>>>>>> })  # 8-9x slower than approach1
>>>>>>>      user  system elapsed
>>>>>>>      9.349   0.280   9.581
>>>>>>> system.time({
>>>>>>> mat2 <- approach2(all.grs[1:100], all.files[1:5])
>>>>>>> })  # 9x slower than previous call, 20x slow than approach1 on same
>>>>>>> input
>>>>>>>      user  system elapsed
>>>>>>>     89.310   0.627  91.044
>>>>>>>           [[alternative HTML version deleted]]
>>>>>>> _______________________________________________
>>>>>>> Bioc-devel at r-project.org mailing list
>>>>>>> https://stat.ethz.ch/mailman/listinfo/bioc-devel
>>>> _______________________________________________
>>>> Bioc-devel at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/bioc-devel

More information about the Bioc-devel mailing list