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
•
written
by
mrodrigues.fernanda •
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.
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!
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
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.
Log in to answer this question.
Could you please post the exact error message?
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.
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