search

dpryan79 / MethylDackel

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

perRead --requireFlags and --ignoreFlags not working properly #117

Closed abcoxyzide closed 3 years ago

abcoxyzide commented 3 years ago

Hi Devon, thanks for developing and maintaining this amazing tool!

I'm playing around with the perRead subcommand, and found out that the --requireFlags and --ignoreFlags do not work as I intended (I expected it to produce the same results as samtools -f / -F). Here's my code:

MethylDackel perRead -@ 2 -q 2 -p 20 --requireFlags 3 --ignoreFlags 3852 ${REF_GENOME} ${IN_FILE} > ${OUT_DIR}/${OUT_FILE}

However, reads with unexpected flags such as 113 and 177 were in the output. If I use samtools view -f3 -F3852, these flags do not show up.

At the moment my work around is to first filter the bam files with samtools and write to disk, then run MethylDackel on those files, since I cannot find a way to pipe stdin to the software.

I'm using version 0.5.2 and htslib-1.11.0

dpryan79 commented 3 years ago

Ah, found the problem. For everything else I have a standard filtering function that gets called. For perRead I used a simplified version of that and only checked if ANY of the flags in the requireFlags were present. I'll fix that in version 0.5.3.