macs3-project/MACS

Q: How to use HMMRATAC to give us a narrow peak file and have signalValue, p and q-value ?

MinaMyr opened this issue · 1 comments

Use case
Hello I want to do peak calling on my ATAC seq data, and for that I use, callpeak and hmmratac.

Describe the problem
My problem here is with hmmratac;

I use this command macs3 hmmratac -f BAM -i RES/atac/try_1/04-NobiasedRegions/ATAC20_D0_rep1_hg38_sort_dedup_biasedRegions.bam -u 15 -l 5 -c 100 -n RES/atac/try_1/06-PeakCalling/hmmratac_narrow/ATAC20_D0_rep1
I read in the doc that we should have a narrow peak file format but I have an *_accessible_regions.gappedPeak format which is not what I want. For the rest of my analysis I need a narrow peak format.

Describe the solution you tried
So I tried to create my narrow peak file from the gappedPeak, by hand, but the trouble I am now having is that for the signal value, p-value and q-value are equals to 0 (I think like this comment said it is intentional #592 (comment) ). Here is what my file look like;

track name="peak" description="peak" type=gappedPeak nextItemButton=on
20	208210	220320	peak_1	0	.	208450	220210	0	3	1240,8820,1550	240,1520,10450	0	0	0
20	220410	267710	peak_2	0	.	220490	267640	0	1	47150	80	0	0	0
20	279100	302420	peak_3	0	.	279110	302390	0	8	890,9080,6700,310,240,3640,1240,660	10,1010,10140,17070,17420,17720,21380,22630	0	0	0
20	302430	314800	peak_4	0	.	302470	314780	0	4	8730,450,2840,130	40,8880,9350,12220	0	0	0
20	314890	342820	peak_5	0	.	314920	342760	0	7	7230,4790,1170,3260,2190,4780,4170	30,7400,12210,13400,16670,18880,23700	0	0	0
20	342840	365480	peak_6	0	.	342880	365440	0	6	5620,10270,120,1640,2240,2450	40,5690,15990,16150,17830,20150	0	0	0
20	368100	402600	peak_7	0	.	368140	402540	0	4	2600,3260,230,28060	40,2650,5950,6380	0	0	0
20	443320	447340	peak_8	0	.	443330	447170	0	2	2860,910	10,2940	0	0	0
20	447420	479680	peak_9	0	.	447440	479630	0	2	15710,16460	20,15750	0	0	0
20	479740	480080	peak_10	0	.	479850	480010	0	1	160	110	0	0	0

I don't really know what to do now, since signal Value p-value and q-value are important for the rest of my analysis.
Thank you in advance !

Hello and thanks for making this resource. I am commenting here to track this since I'm having the same issue: no narrowPeak output from the hmmratac using mapped .bam files with the following code:
macs3 hmmratac \ -i D*rep1_paired.sorted.bam \ *rep2_paired.sorted.bam \ --outdir directory \ -n output_name \ -u 20 \ -l 10 \ -c 2

I previously ran the with the --cutoff-analysis-only flag prior to entering the -u -l -c params. I'm using version 3.0.1 in a conda environment with a SGE scheduler. There are no errors in the output. Here are the last lines of the run output:

INFO @ 14 Jun 2024 17:06:54: [7795 MB] # Write the output... INFO @ 14 Jun 2024 17:09:30: [7795 MB] # Write accessible regions in a gappedPeak file: *accessible_regions.gappedPeak