R encountered fatal error when using processBismarkAln in methylKit
2
0
Entering edit mode
4.7 years ago

Hi

I am attempting to use methylkit to analyse my RRBS data but cannot seem to be able to import my files. I have .bam files generated from bismark. I read I can use function processBismarkAln to read these kind of files into methylkit, but R encounters a fatal error when I try.

Here is my code:

file.list = list("Final_145A.sorted.dedup.bam")
myobj = processBismarkAln(location = file.list, sample.id = list("test145"), assembly = "btaUMD3", save.folder = NULL, save.context = c("CpG"), read.context = "CpG", nolap = FALSE, mincov = 10, minqual = 20, phred64 = FALSE, treatment = c(0))

not sure of what is going on. Do I need to do any extra processing in the bismark .bam file before importing it into methylkit? If yes, what should I do? Specific codes would be appreciated since I am very new to this.

Thank you!!

R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X Yosemite 10.10.5

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] methylKit_1.0.0      devtools_1.12.0      GenomicRanges_1.26.3 GenomeInfoDb_1.10.3 
[5] IRanges_2.8.1        S4Vectors_0.12.1     BiocGenerics_0.20.0 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.9                plyr_1.8.4                 XVector_0.14.0            
 [4] R.methodsS3_1.7.1          bitops_1.0-6               R.utils_2.5.0             
 [7] tools_3.3.3                zlibbioc_1.20.0            mclust_5.2.2              
[10] digest_0.6.12              memoise_1.0.0              tibble_1.2                
[13] gtable_0.2.0               lattice_0.20-34            fastseg_1.20.0            
[16] Matrix_1.2-8               coda_0.19-1                rtracklayer_1.34.2        
[19] withr_1.0.2                stringr_1.2.0              gtools_3.5.0              
[22] Biostrings_2.42.1          grid_3.3.3                 Biobase_2.34.0            
[25] data.table_1.10.4          qvalue_2.6.0               emdbook_1.3.9             
[28] XML_3.98-1.5               BiocParallel_1.8.1         limma_3.30.12             
[31] ggplot2_2.2.1              reshape2_1.4.2             magrittr_1.5              
[34] GenomicAlignments_1.10.0   scales_0.4.1               Rsamtools_1.26.1          
[37] MASS_7.3-45                splines_3.3.3              SummarizedExperiment_1.4.0
[40] assertthat_0.1             bbmle_1.0.18               colorspace_1.3-2          
[43] numDeriv_2016.8-1          stringi_1.1.2              RCurl_1.95-4.8            
[46] lazyeval_0.2.0             munsell_0.4.3              R.oo_1.21.0
methylkit • 2.8k views
ADD COMMENT
0
Entering edit mode

Could you please post the exact error message?

ADD REPLY
0
Entering edit mode

Hi! When R crashed, it just outputted: 'R encountered a fatal errow. Session must be terminated'.

I read somewhere that loading library(Rcpp) before running processBismarkAln could help, but now R is outputting the following:

rsession(53014,0x7fff78526300) malloc: * mach_vm_map(size=8031079741507334144) failed (error code=3) error: can't allocate region ** set a breakpoint in malloc_error_break to debug Error in eval(expr, envir, enclos) : std::bad_alloc Error in eval(expr, envir, enclos) : no methylation tag found.

I am guessing there is something wrong with my .bam file, but not sure how to fix it.

ADD REPLY
0
Entering edit mode

Hello, anyone found a solution for this? I am on windows, and the methylkit crashes when doing processBismarkA : 'R encountered a fatal errow. Session must be terminated'. A colleague uses the same command, the same files, on mac, and it works... Any idea other than switching to MAC (which I am seriously considering...)? Thank you Rita

ADD REPLY
1
Entering edit mode
4.7 years ago

That error message looks like you're not reading in the right bam file - bismark's aligner sets specific flags in the bam file (XM:Z:). If it doesn't encounter these flags it stops (see the code: https://github.com/al2na/methylKit/blob/master/src/methCall.cpp#L560 ) with the error message 'no methylation tag found', which is what you have there.

Judging from your filename, is it possible you ran samtools rmdup on the bismark alignment? This could have removed the bismark tag.

ADD COMMENT
0
Entering edit mode

Hi! Thank you for your input. Because I uses the Ovation RRBS kit from NuGen for library construction, it is recommended to use their customized nudup.py script after alignment to remove duplicates.

Would you have any insights on how to fix the file? I am very new to this. Any help is appreciated.

Thank you!

ADD REPLY
0
Entering edit mode

Looking at the nudup.py script here I cannot see anything where it would remove the tags. Have you tried loading the bam file before you ran nudup.py into methylKit? Does that one work? You may be able to filter the duplicates out in methylkit

ADD REPLY
0
Entering edit mode
4.7 years ago

So that did not work either.Thank you for all your help though, Philipp! The solution I found for the problem was to use the cytosineReport.txt file generated by bismark. However that needs extra processing for loading into methylkit too. If anyone runs through the same problem:

Run the bismark_methylation_extractor; Run the following awk command on the cytosineReport.txt file to put it into the appropriate format for methylkit:

awk '{OFS="\t";if($4+0 > 0 || $5+0 >0 ) print $1,$2,$3,$4/($4+$5),$4+$5;}' cytosineReport.txt > outformethylkit

Then you are able to input these files into methylkit by simply using the read() function. This solution is discussed here: https://groups.google.com/forum/#!topic/methylkit_discussion/0s1DLTWsLyM

I did not seem to find any solution for importing the bismark .bam files into methylkit without R issuing any errors. If anyone has a solution for that, please reply :)

Otherwise, the solution above works.

ADD COMMENT

Login before adding your answer.

Traffic: 1624 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