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