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$start+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]]

