I hope someone can help in this issue.

I have 8 bam files from mm9 alignment, each ~4-5 geg in size. When I run summarizeOverlaps over 3 files, it takes 2-3 hours to finish and it works although my computer almost freezes up. But when I inquire to summarizeOverlaps for the the 8 bam files together, then keep it overnight (as it takes too long to wait), the computer freezes (although it is 16 geg i7 mac, so supposed to be powerful) and the command never results in anything. I even had it run for 30 hours and it looked like it was consuming memory (~600 mega of ram) but still got nothing. I had to reboot the laptop.

I am making my own txdb file from gtf that I used for the alignment to match the naming of the chromosomes. (script is below).

Do you have any tips on how I can get the summzerOverlaps to work on the 8 files to create one se file without freezing up the computer? I have been trying to do that for the past 2 week and always same result.

Any input is appreciated.

here’s the script:

library("DESeq2")
library("GenomicFeatures")
library("Rsamtools")
library("GenomicAlignments")
library("GenomicRanges”)

mm9_from_cluster_gtf_txdb <- makeTranscriptDbFromGFF(file="~/Desktop/genes.gtf", format="gtf”)
head(seqlevels(mm9_from_cluster_gtf_txdb))
saveDb(mm9_from_cluster_gtf_txdb, file="/Path/To/Libraries/TxDB/mm9_from_cluster_Ensembl_txdb.sqlite”)
exonsByGene<-exonsBy(mm9_from_cluster_gtf_txdb,by="gene")
seqinfo(exonsByGene)

fls <- list.files("Path/To/BamFiles", pattern="paired.accepted_hits.bam", full= TRUE)
fls

Experiment <- c(fls[2:8], fls[1])
Experiment

bamLst_experiment <- BamFileList(Experiment, yieldSize=100000)
seqinfo(bamLst_experiment)

se_test_experiment <- summarizeOverlaps(exonsByGene,bamLst_experiment, mode="Union", singleEnd=FALSE, ignore.strand=TRUE, fragments=TRUE) <<<This is the step that freezes the computer when I run the 8 of the files together.

Sessioninfo()

R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base  

other attached packages:
[1] GenomicAlignments_1.2.0 Rsamtools_1.18.1 Biostrings_2.34.0 XVector_0.6.0  
[5] GenomicRanges_1.18.3 GenomeInfoDb_1.2.2 IRanges_2.0.0 S4Vectors_0.4.0  
[9] BiocGenerics_0.12.0 BiocInstaller_1.16.1  

loaded via a namespace (and not attached):
[1] base64enc_0.1-2 BatchJobs_1.5 BBmisc_1.8 BiocParallel_1.0.0 bitops_1.0-6  
[6] brew_1.0-6 checkmate_1.5.0 codetools_0.2-9 DBI_0.3.1 digest_0.6.4  
[11] fail_1.2 foreach_1.4.2 iterators_1.0.7 RSQLite_1.0.0 sendmailR_1.2-1  
[16] stringr_0.6.2 tools_3.1.2 zlibbioc_1.12.0

Hi,

The summarizeOverlaps() method for BamFileList runs the files in parallel with bplapply() from the BiocParallel package. This is described on the man page in the first paragraph of 'Details'.

The default number of cores used in bplapply() is determined by detectCores(), which is generally the max available. While 16GIG was enough to process 3  files concurrently it isn't enough to do all 8 (or the max number of cores you have).

You can control the number of cores used by registering a BPPARAM. The default param on unix and mac is MulticoreParam, on windows it is SnowParam. On my machine you can see the default is a MulticoreParam with 4 workers.

> library(BiocParallel)
> registered()
$MulticoreParam
class: MulticoreParam; bpisup: TRUE; bpworkers: 4; catch.errors: TRUE
setSeed: TRUE; recursive: TRUE; cleanup: TRUE; cleanupSignal: 15;
  verbose: FALSE

$SnowParam
class: SnowParam; bpisup: FALSE; bpworkers: 4; catch.errors: TRUE
cluster spec: 4; type: PSOCK

$BatchJobsParam
class: BatchJobsParam; bpisup: TRUE; bpworkers: NA; catch.errors: TRUE
cleanup: TRUE; stop.on.error: FALSE; progressbar: TRUE

$SerialParam
class: SerialParam; bpisup: TRUE; bpworkers: 1; catch.errors: TRUE

Change the number of workers (and other parameters) by registering a param:

> register(MulticoreParam(workers = 2))
> registered()
$MulticoreParam
class: MulticoreParam; bpisup: TRUE; bpworkers: 2; catch.errors: TRUE
setSeed: TRUE; recursive: TRUE; cleanup: TRUE; cleanupSignal: 15;
  verbose: FALSE

$SnowParam
class: SnowParam; bpisup: FALSE; bpworkers: 4; catch.errors: TRUE
cluster spec: 4; type: PSOCK

$BatchJobsParam
class: BatchJobsParam; bpisup: TRUE; bpworkers: NA; catch.errors: TRUE
cleanup: TRUE; stop.on.error: FALSE; progressbar: TRUE

$SerialParam
class: SerialParam; bpisup: TRUE; bpworkers: 1; catch.errors: TRUE

To run the code in serial,

register(MulticoreParam(workers = 1))

You may have to play around with the number of cores and 'yieldSize' in the BamFileList' to find the optimal combination given the file size, cores and available memory. Let me know how it goes.

Valerie

You might try the featureCounts() function in the Rsubread package. Comparisons published in:

 http://bioinformatics.oxfordjournals.org/content/30/7/923

show that it is considerably faster than compareOverlaps and uses considerably less memory. The efficiency gain is more than an order of magnitude for both time and memory, and is even greater if you have an unusually large number of genomic features. featureCounts() can make use of multiple cores (even without the BiocParallel package).

Hi Valerie-

Thank you for your prompt response. I ran the following 

> library(BiocParallel)
> registered()

OUTPUT:
$MulticoreParam
class: MulticoreParam; bpisup: TRUE; bpworkers: 8; catch.errors: TRUE
setSeed: TRUE; recursive: TRUE; cleanup: TRUE; cleanupSignal: 15; verbose: FALSE

$SnowParam
class: SnowParam; bpisup: FALSE; bpworkers: 8; catch.errors: TRUE
cluster spec: 8; type: PSOCK

$BatchJobsParam
class: BatchJobsParam; bpisup: TRUE; bpworkers: NA; catch.errors: TRUE
cleanup: TRUE; stop.on.error: FALSE; progressbar: TRUE

$SerialParam
class: SerialParam; bpisup: TRUE; bpworkers: 1; catch.errors: TRUE

I will try the script again after using:

register(MulticoreParam(workers = 1))

I will run this overnight now. But, what did you mean by "To run the code in serial .."? Does this mean if i want to run the 8 samples? I am not a computer savvy so i am not sure what most of register() output is indicating initially. After choosing workers=1, i got the following:

$MulticoreParam
class: MulticoreParam; bpisup: TRUE; bpworkers: 1; catch.errors: TRUE
setSeed: TRUE; recursive: TRUE; cleanup: TRUE; cleanupSignal: 15; verbose: FALSE

$SnowParam
class: SnowParam; bpisup: FALSE; bpworkers: 8; catch.errors: TRUE
cluster spec: 8; type: PSOCK

$BatchJobsParam
class: BatchJobsParam; bpisup: TRUE; bpworkers: NA; catch.errors: TRUE
cleanup: TRUE; stop.on.error: FALSE; progressbar: TRUE

$SerialParam
class: SerialParam; bpisup: TRUE; bpworkers: 1; catch.errors: TRUE

 Thank you for your promptness. I've been stuck for almost 2 weeks now at the same step.

By the way, do you recommend that i increase or decrease the yield size. Frankly, the significance of choosing 100,000 was totally arbitrary to me. I did not find a good justification for that.

Thank you Valerie for your feedback. It finally worked. What i did is I used worker=1 and had the run go overnight. It worked without interruptions. I checked it 7 hours later and the results were there (it probably took less time). I maintained the yieldSize as 100000. I guess it is best to run it in serial (not parallel) in case many files would slow it down.

Thank you-

You're welcome. I'm glad it worked for you.

The number of workers to use will depend on file size, the operation and available RAM. You're probably anxious to move to the next step with DESeq2 but I think with some exploration you could use multiple cores. Much of the code in Bioconductor runs in parallel if the resources are available. The general ideas for managing runs are those of yieldSize (or chunking) and choosing the number of workers. You can also use the functions directly from BiocParallel to run your own analysis in parallel. If you have funding available, another resource for large data computing is the Bioconductor AMI:

http://www.bioconductor.org/help/bioconductor-cloud-ami/

Valerie

Thanks Valerie-

I think you make a good point. I will try to run things with 2 workers and yield size 100000 and see how it goes. I will look more into BiocParallel ro run the analysis as well. Thanks a lot!