Open AntonioAhn opened 3 years ago
It seems like the issue is with the third signal track. In state 2, it is 0.187 while in state 1 it is 0.192. I would try even more stringent -u and -l, for instance -u 50 and -l 25. You can also lower the value for --threshold, from 30 to something smaller to increase the peak number - to say 20 or 10. But you will really want to find a set of -u and -l values that correct that issue with the third signal track. That is the likely culprit
Thank you for HMMRATAC. I would love to use this tool but currently, I am having issues with low peak numbers (1262 peaks). This is low in comparison to peaks found with MACS2 (83610 peaks). I have paired-end ATACseq data derived from a cancer cell line.
I ran the tool with:
java -jar HMMRATAC.jar -b ATAC.sorted.bam -i ATAC.sorted.bam.bai -g genome.info -o out_HMM
The model.log shows that the open state (state 2) has the highest emission parameters for all 4 signal components, the nucleosome state (state 1) has the second-highest emission for all signal components and the background state (state 0) has the lowest. I was wondering if my model shows something else wrong that is giving me low peak numbers?
I have also run HMMRATAC with a more stringent setting for creating the model: by increasing -u and -l however this does not change the peak numbers. For example, running with
-u 40 -l 20
gave me 1254 peaks.Would you be able to give some advice on what the issue could be and also how to correct the problem? I have also included a figure showing the fragment size distribution in case this may help.
The model log with default settings shows: