search

dpryan79 / MethylDackel

A (mostly) universal methylation extractor for BS-seq experiments.
MIT License
164 stars 44 forks source link

cytosine_report differences and question #116

Open nchernia opened 3 years ago

nchernia commented 3 years ago

Hello,

I was wondering if there could be an option to stream cytosine_report instead of writing it to an output file. It's pretty big (28G) and I only keep it in order to calculate conversion %.

Also - I've managed to produce a bam file that claims no overlap with any cytosine when it definitely does, and I'm wondering if there's some things getting skipped. I've aligned two ways, once essentially paired-end and once single-end, and the single-end produces overlap statistics but the paired end does not. [There are other slight differences as well, but the reads have the same Phred scores and alignment locations so I'm not sure what's going on.]

This is the call I use: MethylDackel extract -F 1024 -@10 --cytosine_report --CHH --CHG

dpryan79 commented 3 years ago

I suppose that'd doable, maybe a --cytosine_report_stdout option. No clue why the PE data says it has no overlaps, that's quite odd.

nchernia commented 3 years ago

I've sent some test files via email. The "aligning single end separately" version (which I'm calling single) returns a bed file with values when calling extract. The "aligning as paired" one does not - an empty bed file with just the header is what is produced. Interestingly, when I run "perRead" on these same bam files, the problem does not occur and they return the same output, as expected.

nchernia commented 3 years ago

OK, I think this is due to needing the --keepSingleton --keepDiscordant flags.

It would still be nice to have a cytosine report to STD OUT option.

Thanks!