Architecture: considerations on high performance computing with Bioconductor
Overview of performance enhancements
The concept of scalability for software systems and data analyses is complex and precise metrics and criteria are not ready to hand. An interesting overview from the computer science community is that of Weinstock and Goodenough. Scalability of Bioconductor-based analytical workflows can be pursued along various paths, and many useful ideas are presented in work of Mike Lawrence and Martin Morgan, two Bioconductor core members.
This document addresses only the use of parallel computation as directly supported in Bioconductor. The two main approaches are
Shared memory parallelism. The R process is forked an arbitrary number of times with full copies of the memory image and computations proceed independently for each image. Selected results are returned to the master process. This is often effective in multicore environments.
Distributed parallelism. Independent processors, potentially running different operating systems, can run compatible instances of R. Job control can be carried out by R or by a cluster scheduler.
For tasks that are “embarrassingly parallel”, that do not require communication between processes beyond starting, stopping, and returning results, either of these approaches can be used reasonably simply in R.
We can demonstrate the shared memory approach with our laptop.
system.time( lapply(1:8, function(x)Sys.sleep(1) ) )
## user system elapsed ## 0.004 0.001 8.020
##  4
options(mc.cores=4) system.time( mclapply(1:8, function(x)Sys.sleep(1) ) )
## user system elapsed ## 0.007 0.022 2.031
For this meaningless computation, we achieved linear speedup: we cut the time for serial computation by a factor of four.
Detour: BAM files from an RNA-seq experiment
We will be working with a set of BAM files created in the study of Zarnack et al. 2013. The hypothesis of this study is that a class of proteins, the heterogeneous nuclear ribonucleoproteins C1 and C2, referred to as hnRNP C, are responsible for preventing incorporation of Alu elements into mature RNA transcripts. Such incorporation may lead to protein dysfunction that could cause disease; see the references of the Zarnack paper for further background.
The package that we will work with has 8 BAM files with reads aligned to chr14. Four of the files are reads obtained from from HeLa cells (negative control) and four are obtained from HeLa cells in which hnRNP C has been “knocked down” with siRNA.
library(RNAseqData.HNRNPC.bam.chr14) dir(system.file("extdata", package="RNAseqData.HNRNPC.bam.chr14"))
##  "ERR127302_chr14.bam" "ERR127302_chr14.bam.bai" ##  "ERR127303_chr14.bam" "ERR127303_chr14.bam.bai" ##  "ERR127304_chr14.bam" "ERR127304_chr14.bam.bai" ##  "ERR127305_chr14.bam" "ERR127305_chr14.bam.bai" ##  "ERR127306_chr14.bam" "ERR127306_chr14.bam.bai" ##  "ERR127307_chr14.bam" "ERR127307_chr14.bam.bai" ##  "ERR127308_chr14.bam" "ERR127308_chr14.bam.bai" ##  "ERR127309_chr14.bam" "ERR127309_chr14.bam.bai" ##  "tophat2_out"
These are Tabix-indexed BAM files. In the alphabetic ordering, the first four are two pairs of replicates of HeLa cells that have undergone hnRNP C knockdown, and the second four are two pairs of control replicates.
Implicit parallelism through BiocParallel
To foster harmonious development of reliably performant procedures,
Bioconductor introduced the BiocParallel package. Users
benefit from autonomous (but optionally controllable) pursuit
of parallel computing when feasible. Consider the following
example: we will count reads into bins defined by exon addresses
in the HNRNPC example dataset.
The RNAseqData.HNRNPC.bam.chr14 package includes a vector of absolute pathnames of the BAM files, which we assign to
library(RNAseqData.HNRNPC.bam.chr14) fns = RNAseqData.HNRNPC.bam.chr14_BAMFILES length(fns)
##  8
Here we establish the exon bins into which we will count reads.
library(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb = TxDb.Hsapiens.UCSC.hg19.knownGene seqlevels(txdb) = "chr14" ebg = exonsBy(txdb, by="gene")
Now we use the
summarizeOverlaps function from
the GenomicAlignments package to count reads into the exon bins.
We’ll time the counting from a single file, and then time
the counting from all files at once.
library(GenomicAlignments) # summarizeOverlaps uses bplapply to iterate over files s1 = system.time(i1 <- summarizeOverlaps( ebg, fns )) s1
## user system elapsed ## 2.975 0.178 3.178
# show implicit config BiocParallel::bpparam()
## class: MulticoreParam ## bpisup: FALSE; bpnworkers: 4; bptasks: 0; bpjobname: BPJOB ## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE ## bptimeout: 2592000; bpprogressbar: FALSE ## bpRNGseed: ## bplogdir: NA ## bpresultdir: NA ## cluster type: FORK
s2 = system.time(i2 <- summarizeOverlaps( ebg, fns )) s2
## user system elapsed ## 0.271 0.162 10.480
This is not a thorough way of measuring speedup but it
shows reasonable enhancement.
In the second computation, we did approximately 8 times as
much computation, but the clock time elapsed increased only
by a factor of (3.3).
The default behavior of BiocParallel on the MacBook
Air used to produce this document is to pick up the
value of the option
mc.cores and use this as the number
of workers in a
MulticoreParam configuration; if
is NULL, 2 workers are specified.
When BiocParallel is attached,
summarizeOverlaps function will iterate
over the files using
bplapply from the BiocParallel package.
That function will check the R session
for specific parallelization configuration information.
If it finds none, it will check for multiple cores
and make arrangements to use them if present.
The “check” occurs via the function
The default situation on a MacBook Air running MacOSX 10.9.5 with an Intel core i7 processor, (two physical cores with two logical cores each, allowing for four concurrent threads) as follows.
options(mc.cores=NULL) library(BiocParallel) register(MulticoreParam()) bpparam()
## class: MulticoreParam ## bpisup: FALSE; bpnworkers: 2; bptasks: 0; bpjobname: BPJOB ## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE ## bptimeout: 2592000; bpprogressbar: FALSE ## bpRNGseed: ## bplogdir: NA ## bpresultdir: NA ## cluster type: FORK
This identifies an object called
MulticoreParam which is
used to configure the behavior of
bplapply and other utilities
of the BiocParallel package. There are various configuration
object classes that can be used.
‘SnowParam’: distributed memory computing ‘MulticoreParam’: shared memory computing ‘BatchJobsParam’: scheduled cluster computing ‘DoparParam’: foreach computing ‘SerialParam’: non-parallel execution
We need to use
register to determine the type of
concurrent computation that will be performed.
If process size is large, we may want to leave
some cores idle. We can accomplish that by using
Here we’ll check for any advantage of using all logical cores.
library(BiocParallel) register(MulticoreParam(workers=1)) system.time(i3 <- summarizeOverlaps( ebg, fns ))
## user system elapsed ## 15.653 1.286 17.070
register(MulticoreParam(workers=2)) system.time(i3 <- summarizeOverlaps( ebg, fns ))
## user system elapsed ## 0.103 0.029 8.818
register(MulticoreParam(workers=3)) system.time(i3 <- summarizeOverlaps( ebg, fns ))
## user system elapsed ## 0.051 0.023 9.543
register(MulticoreParam(workers=4)) system.time(i3 <- summarizeOverlaps( ebg, fns ))
## user system elapsed ## 0.056 0.026 8.774
all.equal(i3,i2) # check that the results do not change
##  TRUE
Note that in this environment, despite increasing the number of CPUs linearly do not appreciably reduce run time after moving to two cores. This due mainly to communication costs that typically increase with the number of CPUs, and this phenomenon will be environment-specific.
In summary, it is very easy to perform embarrassingly parallel tasks with R, and this carries over to genomic data analysis thanks to BiocParallel. There are some strategic considerations concerning control of memory consumption and communication costs, and full mastery of the topic involves attention to profiling and benchmarking methods that can be addressed in an advanced software development course.