Option: --keepSingleton

[prairie-guy]

First thank you for creating such a useful tool. The speed is awesome! Additionally the flexibility to create C or CpG contexts is very useful.

I am looking to use MD to replace another tool in our pipeline. Before doing so, I am validating the results between the two. At least one issue seems to be associated with how to deal with singletons. For a toy set of data, the only difference between the two tools is a read that has flag 8 set. From your docs, it appears that these are dropped, which is different than the behavior of the other tool.

I had assumed that including the --keepSingleton would solve this, but for me it had no impact. It is still does not seem to report the singleton correctly.

Here is the relevant record from the bam file. The quality scores look good. Perhaps it is being filtered out for another reason.

I would very much appreciate your advice as I would like to validate MD for our pipeline as I have a number of large sequencing files to analyze.

Thanks in advance.

VH01331:238:AAG7FHKM5:1:2610:20182:14953_AACAT  89      17      16497515        60      56M     =       16497515        0       AACCTCCATCTCCCGAAATCAAACGATTCTCCTACCTCAACCTCCCTAATAACTAA      IIIIIIIIIIIIIIIIIIIIIIIIIIIII9IIIIIIIIIIIIIIIIIIIIIIIIII        AS:i:0  NH:i:1  XM:i:0  NM:i:0  MD:Z:16G5G10G5G8G2G2G0G0      YZ:A:-  Yf:i:8  Zf:i:2  XN:i:0  XO:i:0  XG:i:0

[dpryan79]

Your expected behavior is exactly what MethylDackel is doing with --keepSingleton. For reference, the actual code for this is here:

MethylDackel/common.c

Line 429 in 3c77bda

if(!ldata->config->keepSingleton && (b->core.flag & 0x9) == 0x9) continue; //Singleton