vcountPDict Out-of-Memory Errors
0
0
Entering edit mode
@jennylsmith12-16139
Last seen 3.0 years ago

Hello,

I am trying to complete a sequence search using RNA-seq reads with the Biostrings function vcountPDict. However, I am running into out-of-memory errors when running vcountPDict() and I'm looking for some guidance with 1) my implementation may not be optimal or 2) what are typically the largest number of sequences in PDict object one can search for?

Background:

I have 4,722 50bp sequences and thier reverse complement that are saved into a PDict Object (length=9444 sequences). I have on average 5-6 million 75bp RNAseq reads from human Chr 16 that are saved as DNAstringSet object. I am using vcountPDict() function to search the DNAstringSet for the sequences saved in the PDict. But when I run the search with vcountPDict, I get an error "(Error: cannot allocate vector of size 154.7 Gb)" or whatever Gb it was looking for.

I've tried just making the PDict object smaller, but I'm down to 3,000 sequences and I still am getting more than half my searches failing with out-of-memory. So I'd like to hopefully figure out if I'm asking for too many searches and maybe only do 100-200 at a time? Or other suggestions?

CODE:

#create the sequence pattern dictionary
pdict_obj <- sequences %>% 
      pull(SeqExon, name=Name) %>%
      DNAStringSet()
 #Add reverse compliment to search 
 pdict_obj <- c(pdict_obj,  magrittr::set_names(reverseComplement(pdict_obj),
                                      paste0("revComp",names(pdict_obj))))
 pdict_obj <- PDict(pdict_obj)

#Pull out the RNAseq reads for chr16
gr <- GenomicRanges::GRanges(seqnames = Rle("16"),   ranges = IRanges(start=0, end=90354753))
params <- Rsamtools::ScanBamParam(which = gr, what = Rsamtools::scanBamWhat())
RNAseqReads <- Rsamtools::scanBam(bam_file, param = params)[[1]]$seq

#match the sequences
 read.matches <- vcountPDict(pdict = pdict_obj,
                              subject = RNAseqReads,
                              max.mismatch = 0, min.mismatch = 0, 
                              with.indels = FALSE)
  rownames(read.matches) <- names(pdict_obj)

OS Info: I have 80Gb of memory and 4 cores running on Ubuntu version "18.04.4 LTS (Bionic Beaver)".

SessionInfo:

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.4 LTS

Matrix products: default
BLAS/LAPACK: /app/software/OpenBLAS/0.3.7-GCC-8.3.0/lib/libopenblas_haswellp-r0.3.7.so

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                  LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8       LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8      LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

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

other attached packages:
 [1] GGally_2.0.0        gtools_3.8.2        survival_3.2-3      edgeR_3.30.3        limma_3.44.3        RColorBrewer_1.1-2  pheatmap_1.0.12    
 [8] purrr_0.3.4         reshape2_1.4.4      Biostrings_2.56.0   XVector_0.28.0      IRanges_2.22.2      S4Vectors_0.26.1    BiocGenerics_0.34.0
[15] DeGSEA_0.1.0        xlsx_0.6.3          readr_1.3.1         tibble_3.0.1        tidyr_1.1.0         dplyr_1.0.0         gridExtra_2.3      
[22] ggplot2_3.3.2       magrittr_1.5        stringr_1.4.0       knitr_1.29         

loaded via a namespace (and not attached):
 [1] splines_4.0.2          RcppParallel_5.0.2     StanHeaders_2.21.0-5   assertthat_0.2.1       highr_0.8              xlsxjars_0.6.1        
 [7] GenomeInfoDbData_1.2.3 Rsamtools_2.4.0        pillar_1.4.4           lattice_0.20-41        glue_1.4.1             digest_0.6.25         
[13] GenomicRanges_1.40.0   colorspace_1.4-1       Matrix_1.2-18          plyr_1.8.6             pkgconfig_2.0.3        rstan_2.19.3          
[19] zlibbioc_1.34.0        scales_1.1.1           processx_3.4.2         BiocParallel_1.22.0    generics_0.0.2         farver_2.0.3          
[25] ellipsis_0.3.1         withr_2.2.0            cli_2.0.2              crayon_1.3.4           evaluate_0.14          ps_1.3.3              
[31] fansi_0.4.1            pkgbuild_1.0.8         data.table_1.12.8      tools_4.0.2            loo_2.2.0              prettyunits_1.1.1     
[37] hms_0.5.3              lifecycle_0.2.0        matrixStats_0.56.0     munsell_0.5.0          locfit_1.5-9.4         callr_3.4.3           
[43] GenomeInfoDb_1.24.2    compiler_4.0.2         rlang_0.4.6            RCurl_1.98-1.2         grid_4.0.2             rstudioapi_0.11       
[49] bitops_1.0-6           labeling_0.3           gtable_0.3.0           inline_0.3.15          reshape_0.8.8          R6_2.4.1              
[55] rJava_0.9-12           stringi_1.4.6          Rcpp_1.0.5             vctrs_0.3.1            tidyselect_1.1.0       xfun_0.15
software error Biostrings • 636 views
ADD COMMENT

Login before adding your answer.

Traffic: 880 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6