DEXSeq parallelization error
Hi everyone!
I have been attempting to parallelize DEXseq on our SGE cluster using the biocparallel
method described in the latest DEXseq vignette (version 1.16) because it takes forever
to run on one CPU (I have never actually gotten the whole analysis to finish, it gets stuck
at estimateDispersions). I have 26 RNAseq samples with 262,969 exons each.
My Rout report is pasted at the end of this email and my qsub command is as follows:
qsub -pe mpi 5 -l h_vmem=10G -b y -cwd -N dexseq R CMD BATCH dexClust.R
However, once I submit the script to the cluster, it doesn't actually use the 50G of RAM
and it has at its disposal. It uses a maximum mem of 10G and then never actually gets
past the the estimateDispersions step. I am not sure what I am doing wrong. Could you perhaps help? Also, I have tried submitting it to the cluster with one CPU and 20G of RAM,
and it only used 2Gb of RAM....and it has been running for 3 days...
dexClust.Rout:
R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
[Previously saved workspace restored]
> library("DEXSeq", lib.loc="/netapp/home/singh/R/x86_64-unknown-linux-gnu-library/3.2")
Loading required package: BiocParallel
Loading required package: Biobase
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package
arallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':
IQR, mad, xtabs
The following objects are masked from 'package:base':
anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,
do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl,
intersect, is.unsorted, lapply, lengths, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
union, unique, unlist, unsplit
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: IRanges
Loading required package: S4Vectors
Loading required package: stats4
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: DESeq2
Loading required package: SummarizedExperiment
Loading required package: Rcpp
Loading required package: RcppArmadillo
> library("BiocParallel", lib.loc="/netapp/home/singh/R/x86_64-unknown-linux-gnu-library/3.2")
> sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-unknown-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.iso885915 LC_NUMERIC=C
[3] LC_TIME=en_US.iso885915 LC_COLLATE=en_US.iso885915
[5] LC_MONETARY=en_US.iso885915 LC_MESSAGES=en_US.iso885915
[7] LC_PAPER=en_US.iso885915 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.iso885915 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] DEXSeq_1.16.10 DESeq2_1.10.1
[3] RcppArmadillo_0.6.600.4.0 Rcpp_0.12.3
[5] SummarizedExperiment_1.0.2 GenomicRanges_1.22.4
[7] GenomeInfoDb_1.6.3 IRanges_2.4.8
[9] S4Vectors_0.8.11 Biobase_2.30.0
[11] BiocGenerics_0.16.1 BiocParallel_1.4.3
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-2 futile.logger_1.4.1 plyr_1.8.3
[4] XVector_0.10.0 bitops_1.0-6 tools_3.2.1
[7] futile.options_1.0.0 zlibbioc_1.16.0 biomaRt_2.26.1
[10] statmod_1.4.24 rpart_4.1-10 RSQLite_1.0.0
[13] annotate_1.48.0 gtable_0.2.0 lattice_0.20-33
[16] DBI_0.3.1 gridExtra_2.2.1 stringr_1.0.0
[19] hwriter_1.3.2 genefilter_1.52.1 cluster_2.0.3
[22] Biostrings_2.38.4 locfit_1.5-9.1 grid_3.2.1
[25] nnet_7.3-12 AnnotationDbi_1.32.3 XML_3.98-1.4
[28] survival_2.38-3 foreign_0.8-66 latticeExtra_0.6-28
[31] Formula_1.2-1 magrittr_1.5 geneplotter_1.48.0
[34] ggplot2_2.1.0 lambda.r_1.1.7 Rsamtools_1.22.0
[37] Hmisc_3.17-2 scales_0.4.0 splines_3.2.1
[40] colorspace_1.2-6 xtable_1.8-2 stringi_1.0-1
[43] acepack_1.3-3.3 RCurl_1.95-4.8 munsell_0.4.3
> packageVersion("DEXSeq")
[1] '1.16.10'
>
> BPPARAM = MulticoreParam(workers=5)
# dexseq in R
setwd("/netapp/home/singh/pooja_phd/data/pharyngeal_jaw_transcriptomes/raw_demultiplexed_bam/paper1/map2Onil/all_bam/alternative_splicing/")
>
>
> # reads files
> countFiles = list.files(".", pattern="DEXseqcount.txt$", full.names=TRUE)
> countFiles
[1] "./31246_mapped.sorted.PE.dedup.DEXseqcount.txt"
[2] "./31247_mapped.sorted.PE.dedup.DEXseqcount.txt"
[3] "./31248_mapped.sorted.PE.dedup.DEXseqcount.txt"
[4] "./31249_mapped.sorted.PE.dedup.DEXseqcount.txt"
[5] "./31250_mapped.sorted.PE.dedup.DEXseqcount.txt"
[6] "./31251_mapped.sorted.PE.dedup.DEXseqcount.txt"
[7] "./31252_mapped.sorted.PE.dedup.DEXseqcount.txt"
[8] "./31253_mapped.sorted.PE.dedup.DEXseqcount.txt"
[9] "./Index04_mapped.sorted.PE.dedup.DEXseqcount.txt"
[10] "./Index05_mapped.sorted.PE.dedup.DEXseqcount.txt"
[11] "./Index06_mapped.sorted.PE.dedup.DEXseqcount.txt"
[12] "./Index07_mapped.sorted.PE.dedup.DEXseqcount.txt"
[13] "./Index08_mapped.sorted.PE.dedup.DEXseqcount.txt"
[14] "./Index09_mapped.sorted.PE.dedup.DEXseqcount.txt"
[15] "./Index10_mapped.sorted.PE.dedup.DEXseqcount.txt"
[16] "./Index11_mapped.sorted.PE.dedup.DEXseqcount.txt"
[17] "./Index12_mapped.sorted.PE.dedup.DEXseqcount.txt"
[18] "./Index13_mapped.sorted.PE.dedup.DEXseqcount.txt"
[19] "./Index14_mapped.sorted.PE.dedup.DEXseqcount.txt"
[20] "./Index15_mapped.sorted.PE.dedup.DEXseqcount.txt"
[21] "./Index16_mapped.sorted.PE.dedup.DEXseqcount.txt"
[22] "./Index18_mapped.sorted.PE.dedup.DEXseqcount.txt"
[23] "./Index19_mapped.sorted.PE.dedup.DEXseqcount.txt"
[24] "./Index20_mapped.sorted.PE.dedup.DEXseqcount.txt"
[25] "./Index21_mapped.sorted.PE.dedup.DEXseqcount.txt"
>
> #read gff for DEXseq
> flattenedFile = list.files("/netapp/home/singh/pooja_phd/data/cichlid_reference_genomes/new_Oreochromis/gtf_brawand/", pattern="seq.gff$", full.names=TRUE)
> flattenedFile
[1] "/netapp/home/singh/pooja_phd/data/cichlid_reference_genomes/new_Oreochromis/gtf_brawand//Oreochromis_niloticus.BROADON2.DEXseq.gff"
>
> # import sample table
> sampleTable <- read.table("paper1_speciesinfo.txt", header=T, row.names=1)
>
> #start analysis
> dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable, design= ~ sample + exon + condition:exon, flattenedfile=flattenedFile)
converting counts to integer mode
>
>
> # normlisation
>
> dxd = estimateSizeFactors(dxd)
>
> dxd = estimateDispersions(dxd, BPPARAM = BPPARAM)
I would really appreciate any suggestions anyone might have!
I have been attempting to parallelize DEXseq on our SGE cluster using the biocparallel
method described in the latest DEXseq vignette (version 1.16) because it takes forever
to run on one CPU (I have never actually gotten the whole analysis to finish, it gets stuck
at estimateDispersions). I have 26 RNAseq samples with 262,969 exons each.
My Rout report is pasted at the end of this email and my qsub command is as follows:
qsub -pe mpi 5 -l h_vmem=10G -b y -cwd -N dexseq R CMD BATCH dexClust.R
However, once I submit the script to the cluster, it doesn't actually use the 50G of RAM
and it has at its disposal. It uses a maximum mem of 10G and then never actually gets
past the the estimateDispersions step. I am not sure what I am doing wrong. Could you perhaps help? Also, I have tried submitting it to the cluster with one CPU and 20G of RAM,
and it only used 2Gb of RAM....and it has been running for 3 days...
dexClust.Rout:
R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
[Previously saved workspace restored]
> library("DEXSeq", lib.loc="/netapp/home/singh/R/x86_64-unknown-linux-gnu-library/3.2")
Loading required package: BiocParallel
Loading required package: Biobase
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package
arallel':clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':
IQR, mad, xtabs
The following objects are masked from 'package:base':
anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,
do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl,
intersect, is.unsorted, lapply, lengths, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
union, unique, unlist, unsplit
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: IRanges
Loading required package: S4Vectors
Loading required package: stats4
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: DESeq2
Loading required package: SummarizedExperiment
Loading required package: Rcpp
Loading required package: RcppArmadillo
> library("BiocParallel", lib.loc="/netapp/home/singh/R/x86_64-unknown-linux-gnu-library/3.2")
> sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-unknown-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.iso885915 LC_NUMERIC=C
[3] LC_TIME=en_US.iso885915 LC_COLLATE=en_US.iso885915
[5] LC_MONETARY=en_US.iso885915 LC_MESSAGES=en_US.iso885915
[7] LC_PAPER=en_US.iso885915 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.iso885915 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] DEXSeq_1.16.10 DESeq2_1.10.1
[3] RcppArmadillo_0.6.600.4.0 Rcpp_0.12.3
[5] SummarizedExperiment_1.0.2 GenomicRanges_1.22.4
[7] GenomeInfoDb_1.6.3 IRanges_2.4.8
[9] S4Vectors_0.8.11 Biobase_2.30.0
[11] BiocGenerics_0.16.1 BiocParallel_1.4.3
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-2 futile.logger_1.4.1 plyr_1.8.3
[4] XVector_0.10.0 bitops_1.0-6 tools_3.2.1
[7] futile.options_1.0.0 zlibbioc_1.16.0 biomaRt_2.26.1
[10] statmod_1.4.24 rpart_4.1-10 RSQLite_1.0.0
[13] annotate_1.48.0 gtable_0.2.0 lattice_0.20-33
[16] DBI_0.3.1 gridExtra_2.2.1 stringr_1.0.0
[19] hwriter_1.3.2 genefilter_1.52.1 cluster_2.0.3
[22] Biostrings_2.38.4 locfit_1.5-9.1 grid_3.2.1
[25] nnet_7.3-12 AnnotationDbi_1.32.3 XML_3.98-1.4
[28] survival_2.38-3 foreign_0.8-66 latticeExtra_0.6-28
[31] Formula_1.2-1 magrittr_1.5 geneplotter_1.48.0
[34] ggplot2_2.1.0 lambda.r_1.1.7 Rsamtools_1.22.0
[37] Hmisc_3.17-2 scales_0.4.0 splines_3.2.1
[40] colorspace_1.2-6 xtable_1.8-2 stringi_1.0-1
[43] acepack_1.3-3.3 RCurl_1.95-4.8 munsell_0.4.3
> packageVersion("DEXSeq")
[1] '1.16.10'
>
> BPPARAM = MulticoreParam(workers=5)
# dexseq in R
setwd("/netapp/home/singh/pooja_phd/data/pharyngeal_jaw_transcriptomes/raw_demultiplexed_bam/paper1/map2Onil/all_bam/alternative_splicing/")
>
>
> # reads files
> countFiles = list.files(".", pattern="DEXseqcount.txt$", full.names=TRUE)
> countFiles
[1] "./31246_mapped.sorted.PE.dedup.DEXseqcount.txt"
[2] "./31247_mapped.sorted.PE.dedup.DEXseqcount.txt"
[3] "./31248_mapped.sorted.PE.dedup.DEXseqcount.txt"
[4] "./31249_mapped.sorted.PE.dedup.DEXseqcount.txt"
[5] "./31250_mapped.sorted.PE.dedup.DEXseqcount.txt"
[6] "./31251_mapped.sorted.PE.dedup.DEXseqcount.txt"
[7] "./31252_mapped.sorted.PE.dedup.DEXseqcount.txt"
[8] "./31253_mapped.sorted.PE.dedup.DEXseqcount.txt"
[9] "./Index04_mapped.sorted.PE.dedup.DEXseqcount.txt"
[10] "./Index05_mapped.sorted.PE.dedup.DEXseqcount.txt"
[11] "./Index06_mapped.sorted.PE.dedup.DEXseqcount.txt"
[12] "./Index07_mapped.sorted.PE.dedup.DEXseqcount.txt"
[13] "./Index08_mapped.sorted.PE.dedup.DEXseqcount.txt"
[14] "./Index09_mapped.sorted.PE.dedup.DEXseqcount.txt"
[15] "./Index10_mapped.sorted.PE.dedup.DEXseqcount.txt"
[16] "./Index11_mapped.sorted.PE.dedup.DEXseqcount.txt"
[17] "./Index12_mapped.sorted.PE.dedup.DEXseqcount.txt"
[18] "./Index13_mapped.sorted.PE.dedup.DEXseqcount.txt"
[19] "./Index14_mapped.sorted.PE.dedup.DEXseqcount.txt"
[20] "./Index15_mapped.sorted.PE.dedup.DEXseqcount.txt"
[21] "./Index16_mapped.sorted.PE.dedup.DEXseqcount.txt"
[22] "./Index18_mapped.sorted.PE.dedup.DEXseqcount.txt"
[23] "./Index19_mapped.sorted.PE.dedup.DEXseqcount.txt"
[24] "./Index20_mapped.sorted.PE.dedup.DEXseqcount.txt"
[25] "./Index21_mapped.sorted.PE.dedup.DEXseqcount.txt"
>
> #read gff for DEXseq
> flattenedFile = list.files("/netapp/home/singh/pooja_phd/data/cichlid_reference_genomes/new_Oreochromis/gtf_brawand/", pattern="seq.gff$", full.names=TRUE)
> flattenedFile
[1] "/netapp/home/singh/pooja_phd/data/cichlid_reference_genomes/new_Oreochromis/gtf_brawand//Oreochromis_niloticus.BROADON2.DEXseq.gff"
>
> # import sample table
> sampleTable <- read.table("paper1_speciesinfo.txt", header=T, row.names=1)
>
> #start analysis
> dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable, design= ~ sample + exon + condition:exon, flattenedfile=flattenedFile)
converting counts to integer mode
>
>
> # normlisation
>
> dxd = estimateSizeFactors(dxd)
>
> dxd = estimateDispersions(dxd, BPPARAM = BPPARAM)
I would really appreciate any suggestions anyone might have!