Thanks Mike and Martin for your advice. >From within the script, the .Rout didn't report the file name in the second case. It turned out to be two different errors in two different files, which I identified using the trace you suggested. Other files don't give that problem. Below my tests and the output: I cut repeated text, but it's still quite long. My apologies for that. I run Picard ValidateSam on several files. there are a few errors of this kind ERROR: Record 57116264, Read name HS16_6904:5:1203:11452:139702#8, MAPQ should be 0 for unmapped read. with the only difference that the file kc5165129.bam generating the "Realloc" error has a bit more of those, compared to the one generating a segmentation fault. these few errors shouldn't be the problem however, because they are present also in the files which run smoothly. thanks, Francesco ###################### # first test ######################## trace(scanBam, quote({ print(file); print(param) })) some.reads <- scanBam("kc5165129.bam", param=ScanBamParam(what=c("pos","mapq"), which=GRanges("1",IRanges(1,1000000)))) some.reads <- scanBam("kc5165122.bam", param=ScanBamParam(what=c("pos","mapq"), which=GRanges("1",IRanges(1,1000000)))) open(BamFile("kc5165129.bam", "kc5165129.bam"), "rb") ################# # full test ################## trace(scanBam, quote({ print(file); print(param) })) target.df <- read.delim("/share/apps/genomics/reference/capture/refseq _exonshg19_exons.bed.dict.merged.bed", header=FALSE,col.names=c("seqname","start","end","exon")) target <- GRanges(seqname=target.df$seqname,IRanges(start=target.df$st art+1,end=target.df$end)) target <- reduce(target) target.sub <- subdivideGRanges(target) test.count <- countBamInGRanges("kc5165129.bam",target.sub) Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_3/kc5165129.bam index: /mypath/ICH2/k/bam_batch_3/kc5165129.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq # #### [ all similar entries - cut to simplify report ] # Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_3/kc5165129.bam index: /mypath/ICH2/k/bam_batch_3/kc5165129.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_3/kc5165129.bam index: /mypath/ICH2/k/bam_batch_3/kc5165129.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq Error in value[[3L]](cond) : Realloc could not re-allocate memory (20971520 bytes) file: /mypath/ICH2/k/bam_batch_3/kc5165129.bam index: /mypath/ICH2/k/bam_batch_3/kc5165129.bam #################### # test another file #################### test.count <- countBamInGRanges("kc5165122.bam",target.sub) Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_3/kc5165122.bam index: /mypath/ICH2/k/bam_batch_3/kc5165122.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq *** caught segfault *** address (nil), cause 'memory not mapped' Traceback: 1: .Call(.NAME, ..., PACKAGE = PACKAGE) 2: .Call2(func, .extptr(file), space, flag, simpleCigar, ..., PACKAGE = "Rsamtools") 3: doTryCatch(return(expr), name, parentenv, handler) 4: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 5: tryCatchList(expr, classes, parentenv, handlers) 6: tryCatch({ .Call2(func, .extptr(file), space, flag, simpleCigar, ..., PACKAGE = "Rsamtools")}, error = function(err) { stop(conditionMessage(err), "\n file: ", path(file), "\n index: ", index(file))}) 7: .io_bam(.scan_bamfile, file, param = param, reverseComplement, tmpl) 8: scanBam(bam, param = param) 9: scanBam(bam, param = param) 10: eval(expr, envir, enclos) 11: eval(call, sys.frame(sys.parent())) 12: callGeneric(bam, ..., param = param) 13: scanBam(bam.file, param = ScanBamParam(what = scan.what, which = range(granges.subset))) 14: scanBam(bam.file, param = ScanBamParam(what = scan.what, which = range(granges.subset))) 15: countBamInGRanges("kc5165122.bam", target.sub) Possible actions: 1: abort (with core dump, if enabled) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection: 1 aborting ... Segmentation fault ####################### # test a third file ####################### test.count <- countBamInGRanges("kc5165146.bam",target.sub) Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_4/kc5165146.bam index: /mypath/ICH2/k/bam_batch_4/kc5165146.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq # #### [ all similar entries - cut to simplify report ] # Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_4/kc5165146.bam index: /mypath/ICH2/k/bam_batch_4/kc5165146.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq Tracing scanBam(bam, param = param) on entry class: BamFile path: /mypath/ICH2/k/bam_batch_4/kc5165146.bam index: /mypath/ICH2/k/bam_batch_4/kc5165146.bam isOpen: TRUE class: ScanBamParam bamFlag (NA unless specified): bamSimpleCigar: FALSE bamReverseComplement: FALSE bamTag: bamWhich: 24 elements bamWhat: pos, mapq On 18 Sep 2012, at 17:24, Martin Morgan <mtmorgan@fhcrc.org<mailto:mtmorgan@fhcrc.org>> wrote: On 09/18/2012 07:43 AM, Lescai, Francesco wrote: Hi there, I'm encountering a memory error with exomeCopy which might be file- specific (i.e. it's not happening with any of my other batches). I'm running it on a cluster. Initially I had three batches pooled together and I had this error message in my .Rout file ---------------- for (f1 in bamfiles){ + sample<-gsub('.bam','',f1) + exome.samples<-c(exome.samples,sample) + rdata[[sample]] <- countBamInGRanges(f1,target.sub) + } Error in value[[3L]](cond) : Realloc could not re-allocate memory (20971520 bytes) file: /SAN/biomed/biomed14/ICH2/k/bam_batch_3/K_5165129.bam index: /SAN/biomed/biomed14/ICH2/k/bam_batch_3/K_5165129.bam Calls: countBamInGRanges ... tryCatch -> tryCatchList -> tryCatchOne -> <anonymous> Execution halted ------------------ The index doesn't seem to point to a .bam.bai file, but I checked of course and the .bai is there, same way as for the other bams I'm counting in. The index uses the samtools naming convention, so '.bam' without the explicit '.bai'. The Realloc error could indicate out-of-memory, or possibly a bug (but I would have expected an error about not able to allocate 0 bytes) that is fixed in the devel (soon to be release) version of Rsamtools http://bioconductor.org/developers/useDevel/ Then I separated the batches, and run again. Two of three didn't display any error, while the third one (the one that .bam file belongs to) had this (slightly different) error ------------------------- for (f1 in bamfiles){ + sample<-gsub('.bam','',f1) + exome.samples<-c(exome.samples,sample) + rdata[[sample]] <- countBamInGRanges(f1,target.sub) + } *** caught segfault *** address 0x4, cause 'memory not mapped' Traceback: 1: .Call(.bamfile_open, path(con), index(con), "rb") 2: open.BamFile(BamFile(file, index), "rb") 3: open(BamFile(file, index), "rb") 4: scanBam(bam.file, param = ScanBamParam(what = scan.what, which = range(granges.subset))) 5: scanBam(bam.file, param = ScanBamParam(what = scan.what, which = range(granges.subset))) 6: countBamInGRanges(f1, target.sub) aborting ... ------------------------- Any idea where I could look better, in order to track this error? but this looks like a different error, unrelated to memory allocation. Mike's suggestion to narrow it to invoking essentially step 3 in the traceback -- open(BamFile(file, index), "rb") -- is the best way forward. To get there, you might say trace(scanBam, quote({ print(file); print(param) })) to figure out where in countBamInGRanges() things are going wrong. Martin thanks for your help, Francesco sessionInfo() R version 2.15.0 (2012-03-30) 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=C LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.iso885915 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] exomeCopy_1.2.0 Rsamtools_1.8.5 Biostrings_2.24.1 [4] GenomicRanges_1.8.7 IRanges_1.14.4 BiocGenerics_0.2.0 loaded via a namespace (and not attached): [1] bitops_1.0-4.1 stats4_2.15.0 zlibbioc_1.2.0 [[alternative HTML version deleted]] _______________________________________________ Bioconductor mailing list Bioconductor@r-project.org<mailto:bioconductor@r-project.org> https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor -- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793 ---------------------------------------------------------------------- ----------- Francesco Lescai, PhD, EDBT Senior Research Associate in Genome Analysis University College London Faculty of Population Health Sciences Dept. Genes, Development & Disease ICH - Molecular Medicine Unit, GOSgene team 30 Guilford Street WC1N 1EH London UK email: f.lescai@ucl.ac.uk<mailto:f.lescai@ucl.ac.uk> phone: +44.(0)207.905.2274 [ext: 2274] ---------------------------------------------------------------------- ---------- [[alternative HTML version deleted]]