[BioC] Caught segfault when running MEDIPS package

[Asta Laiho [guest]]

Hi,

I'm running MEDIPS package to analyze methylation sequencing data. I have many samples (from the same source and preprocessed the same way) and some of them are performing without problems but for some there is an error with segfault during the methylProfiling step. For some of the problematic samples it helps when I manually reload the saved normalized data object and rerun the code interactively but for one of the samples I still get the same error. 

The code, error and session info are given below. I would be very grateful for any help!


genome <- c("BSgenome.Hsapiens.UCSC.hg19")
pgenome <- "hg19"
library(BSgenome.Hsapiens.UCSC.hg19)

coverage.resolution <- 50
smoothing.extension <- 400

dat <- MEDIPS.readAlignedSequences(BSgenome=genome, file=paste(aDir, "/MEDIPS_input/",samples[i], "_MEDIPS-input.tsv", sep=""))
dat <- MEDIPS.genomeVector(data=dat, bin_size=as.numeric(coverage.resolution), extend=as.numeric(smoothing.extension))
dat <- MEDIPS.getPositions(data=dat, pattern="CG")
dat <- MEDIPS.couplingVector(data=dat, fragmentLength=as.numeric(fragment.length), func="count")
dat <- MEDIPS.normalize(data=dat)


roisFile <- paste(aDir, "/promoter_rois.xls", sep="")
if(!file.exists(roisFile)){
 library(rtracklayer)
 session <- browserSession()
 genome(session) <- pgenome
 query <- ucscTableQuery(session, "refGene")
 refGene <- getTable(query)
 sta <- refGene$txStart
 sta[refGene$strand=="+"] <- sta[refGene$strand=="+"]-1000
 sta[refGene$strand=="-"] <- sta[refGene$strand=="-"]-500
 sto <- refGene$txStart
 sto[refGene$strand=="+"] <- sta[refGene$strand=="+"]+500
 sto[refGene$strand=="-"] <- sta[refGene$strand=="-"]+1000
 rois <- unique(data.frame(refGene$chrom, sta, sto, make.unique(as.character(refGene$name))))
 write.table(rois, file=roisFile, sep="\t", quote=F, col.names=F, row.names=F)
rois <- read.table(roisFile, sep="\t", as.is=T)
colnames(rois) <- c("refGene.chrom", "sta", "sto", "refGene.name")

frames <- MEDIPS.methylProfiling(data1=dat, ROI_file=roisFile, math=mean, select=2)

Preprocessing...
Reading ROIs...
Extract data according to given ROI...
Methylation profile will be calculated on the ROI data set
Analysed 39041 / 41310 
 *** caught segfault ***
address 0x2aca1311f900, cause 'memory not mapped'

Traceback:
 1: .Call("roiprofile", input, as.numeric(select), as.matrix(ROI2),     as.integer(chr_binposition), data1, data2, environment(wilcox.test),     wilcox.test, environment(var), var, environment(math), math,     t.test, environment(t.test), as.numeric(factor(chr_names(data1))))
 2: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
 3: suppressWarnings(.Call("roiprofile", input, as.numeric(select),     as.matrix(ROI2), as.integer(chr_binposition), data1, data2,     environment(wilcox.test), wilcox.test, environment(var),     var, environment(math), math, t.test, environment(t.test),     as.numeric(factor(chr_names(data1)))))
 4: MEDIPS.methylProfiling(data1 = dat, ROI_file = roisFile, math = mean,     select = 2)

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:


 -- output of sessionInfo(): 

R version 2.14.0 (2011-10-31)
Platform: x86_64-redhat-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] 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   
 [7] LC_PAPER=C                 LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] MEDIPS_1.4.0                       BSgenome.Hsapiens.UCSC.hg19_1.3.17
[3] BSgenome_1.22.0                    Biostrings_2.22.0                 
[5] GenomicRanges_1.6.4                IRanges_1.12.5                    

loaded via a namespace (and not attached):
[1] gtools_2.6.2


--
Sent via the guest posting facility at bioconductor.org.