Getting GC context methylation report

[erlitzkin]

Hello,

I am running into an issue when I try running coverage2cytosine with the --gc/--gc_context option.
I run the following command:
coverage2cytosine --gc_context --genome_folder <path> -o <output> [input.cov.gz]

My input file is the *cov.gz file generated by first running bismark_methylation_extractor with the CX context option:
bismark_methylation_extractor --report --buffer_size 10G --gzip --CX --bedGraph --counts --output <output> [input.bam]

When I run the coverage2cytosine GC command, this is what happens:
The first part, "Methylation information will now be written into a genome-wide cytosine report," works as expected. However, the second part, "Methylation information for GC context will now be written to a GpC-context report" fails and returns the following error:
Illegal division by zero at /path/to/opt/anaconda3/bin/coverage2cytosine line 984, <IN> line 48318.
(In the above, I changed my directory names to "/path/to" for privacy reasons).

I would greatly appreciate any help in resolving this, as I have been unable to find a solution on my own so far.

[FelixKrueger]

Hmm, I would assume that your command is correct, I am not sure if you also need to add --CX for the coverage2cytosine command?

Other things:

• which version are you using?
• could you check the input file for signs of truncation? Something like this might help: gunzip -c input.cov.gz | head -48323 | tail

[erlitzkin]

Hi Felix,

Thanks so much for your speedy reply!

I am using v0.23.1.

This is what I get when I run the gunzip command you gave in your reply:

When you say to add --CX to the command, to you mean like this?
coverage2cytosine --CX --gc_context --genome_folder <path> -o <output> [input.cov.gz]
When I run the above command, I get the same error.

[FelixKrueger]

Thanks for spotting this, this bug was triggered when the coverage was set to 0. I have now fixed it in the dev version (here), can you clone the dev branch and try again?

[erlitzkin]

Hi Felix,

Excellent news: it seems to be working now! Apologies that it took some time for me to test it out -- this was my first time cloning a dev branch or running any program from within it, so it took me a minute to figure out (I am learning!). Thank you so much for your help and for debugging.

Best,
Noa

[FelixKrueger]

Excellent, glad it's working well now!